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UX.tiY  The  derivative  of  the  local  u  velocity  with  respect  to 

x  and  y  evaluated  in  the  computational  plane 


v 


Variable  which  denotes  the  local  y  component  of  the  flow 
field  velocity 
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LIST  OF  SYMBOLS  (Continued) 


SYMBOL 

v*  Current  value  of  the  v  velocity  component  obtained  from 

the  finite  difference  equation 

v  Normal  field  velocity  component  relative  to  the  surface 

normal  direction 

V  Control  volume 

vk’^k  The  ^  velocity  component  in  indicial  notation 

VC  Quantity  which  is  proportional  to  the  component  of  the 

local  fluid  velocity  in  the  n  direction 

2 

VV  Combination  of  grid  transformation  coefficients  -o/2J  . 

VX,VY  The  derivative  of  the  local  v  velocity  component  with 

respect  to  x  and  y  evaluated  in  the  computational  plane 

w  Function  which  gives  the  maximum  velocity  defect  in  the 

wake  at  each  location  downstream 

W  Nondimensional  complex  valued  velocity  u-iv 

W  Complex  valued  velocity  (u-iv) 

WL,  Non-interacting  length  in  the  near  wake  measured  from 

the  trailing  edge  of  the  airfoil 

MLg  Turbulent  interaction  length  in  the  near  wake  measured 

from  the  trailing  edge  of  the  airfoil 

WLav  Average  value  of  WL-j  and  WL,, 

x  Coordinate  in  the  physical  plane  defining  the  airfoil 

chord  axis;  and  a  local  tangential  coordinate  in  the 
boundary-layer  turbulence  model 

x^  The  kth  coordinate  in  indicial  notation 

x.  Turbulent  transition  length  measured  along  the  body 

surface 

x  The  x  coordinate  of  a  location  in  the  physical  plane 

"  about  which  a  moment  is  computed;  and  a  new  x  coordinate 

location  obtained  by  the  interpolation  of  adjacent  £ 
grid  lines 
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LIST  OF  SYMBOLS  (Continued) 


SYMBOL 

xR  The  x  coordinate  which  locates  the  position  of  the  shear 

layer  reattachment  which  forms  the  bubble  in  the  numerical 
solution 

Xx  The  x  coordinate  which  locates  the  start  of  transition  to 

turbulence  on  the  upper  airfoil  surface  in  the  numerical 
solution 

x-j-  The  x  coordinate  which  locates  where  fully  turbulent  flow 

first  occurs  downstream  of  the  leading  edge  on  the  upper 
airfoil  surface  in  the  numerical  solution 

XETA  Modified  grid  transformation  metric  x^/Zj 

XXI  Modified  grid  transformation  metric  x,/2j 

y  Coordinate  in  the  physical  plane  normal  to  the  airfoil 

chord;  a  local  normal  coordinate  in  the  boundary-layer 
turbulence  model 

yk  Computed  location  of  the  kth  point  within  the  boundary- 

'  layer  for  equal  velocity  increments  in  the  coordinate 

attraction  method 

y  The  y  coordinate  of  a  location  in  the  physical  plane  about 

p  which  a  moment  is  computed;  a  new  y  coordinate  location 

obtained  by  the  interpolation  of  adjacent  5  grid  lines 

y.(x)  Airfoil  thickness  function  for  NACA  four  digit  airfoil 

sections 

YETA  Modified  grid  transformation  metric  y^/2j 

YXI  Modified  grid  transformation  metric  y^./2j 

z  Input  axis  variable  for  calculating  x  coordinates  of  the 

final  distributed  surface  points 

Z  Nondimensional  complex  variable  x+iy  in  the  far  field 

boundary  condition  model 

Z  Complex  variable  c(x+iy)  in  the  far  field  boundary  con¬ 

dition  model 

a  Geometric  angle  of  attack;  and  a  grid  transformation 

coefficient  x2  +  y2. 

n  n 
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LIST  OF  SYMBOLS  (Continued) 

SYMBOL 

Geometric  angle  of  attack 

a-j  The  value  f 1 '  (0)  in  the  Blasius'  series. 

8  Grid  transformation  coefficient  x_x  + 

i  n  ■'C  n 

Y  Grid  transformation  coefficient  x*  +  y|;  and  the 

intermittency  factor  defined  in  Equation  63 

r  Clockwise  circulation  around  a  body 

r(x)  Turbulent  transition  factor  function 

6  Local  boundary-layer  thickness 

6  Local  wake  half  width 

w 

6*  Local  displacement  thickness  defined  by  Equation  64 

Kronecker  Delta 

J  ’ 

Ay^  Input  displacement  in  the  last  two  y  positions  at  the  outer 

boundary  in  the  coordinate  line  attraction  method 

e  Dissipation  rate  of  turbulent  energy 

e .  Inner  layer  eddy  viscosity 

Turbulent  eddy  viscosity 

eQ  Outer  layer  eddy  viscosity 

ew  Wake  value  of  the  eddy  viscosity 

e.  Inner  layer  eddy  viscosity  for  comparison  in  limiting 

technique 

eQW  Far  wake  eddy  viscosity  value 

n  Coordinate  in  the  computational  plane  which  is  used 

to  define  the  body  surface 
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LIST  OF  SYMBOLS  (Continued) 


SYMBOL 

n Summation  parameter  defined  as  (k-l)An 

n  A  selected  contour  integral  path  around  the  body 

e  Local  momentum  thickness;  and  the  angular  coordinate 

in  complex  polar  variables 

0(5)  Function  which  specifies  the  angle  that  the  incoming 

£  lines  make  with  the  body  surface  coordinate 

A  Transition  distance  defined  as  the  distance  required  to 

vary  r  (x)  from  1/4  to  3/4 

Ad  Turbulent  reattachment  criteria  defined  as  (0/u  )  du  /ds 

K  G  6 

n  Dynamic  viscosity  coefficient 

lip  Potential  doublet  strength 

v  Kinematic  viscosity  u/p 


£ 


IT 

P 

a 


aB’°n 

T 

Tij 
Tti  j 
THj 


Coordinate  in  the  computational  plane  which  emanates  from 
the  body  surface  grid  points  and  is  orthogonal  to  n  in 
the  computational  plane 

Ratio  of  circumference  to  diameter  of  a  circle 
Fluid  density 

Similarity  parameter  in  Blaslus1  series  of  Equation  16; 
and  a  grid  transformation  coefficient  defined  by  Equation 
A- 7 


Inner  and  outer  control  volume  surfaces 

Grid  transformation  coefficient  defined  by  Equation  A-8 
Element  of  the  viscous  stress  tensor 

Element  of  the  turbulent  stress  tensor 

Element  of  the  total  stress  tensor, 

Transonic  small  disturbance  potential  doublet  function 
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LIST  OF  SYMBOLS  (Continued) 

SYMBOL 

$  Transonic  small  disturbance  potential  vortex  function 

$  Acceleration  parameter  for  pressure  iteration  given  by 

2 

Equation  40  defined  as  2^^  /  ((a  +  y)At) 
if>  Inviscid  stream  function 

u  Magnitude  of  the  local  flow  field  vorticity,  curl  v_; 

and  a  local  acceleration  parameter 

id  Acceleration  parameter  for  the  field  pressure  finite 

^  difference  equations 

ojpb  Acceleration  parameter  for  the  surface  pressure  iteration 

id  Acceleration  parameter  for  the  u  and  v  velocity  component 

finite  difference  equations 


Denotes  a  location  on  the  body  surface 


e  Denotes  a  location  at  the  edge  of  the  boundary-layer  or 

wake 


f  Denotes  a  location  on  the  computational  far  field  boundary 

i  Denotes  the  ?  coordinate  location  in  the  computational 

plane  for  finite  difference  equations;  and  an  indicial 
notation  index  in  differential  equations 

j  Denotes  the  n  coordinate  location  in  the  computational 

plane  for  finite  difference  equations;  and  an  indicial 
notation  index  in  differential  equations 

k  A  component  index  in  indicial  notation 

max, min  Denote  the  maximum  and  minimum  values  of  a  variable 

S  Denotes  a  value  at  the  laminar  separation  point  which 

defines  the  beginning  of  the  bubble 

w  Denotes  a  value  in  the  wake  region 

x,xx  Denote  differentiation  with  respect  to  x 
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SYMBOL 

y.yy 

n  »nn 


LIST  OF  SYMBOLS  (Concluded) 

Denote  differentiation  with  respect  to  y 

Denote  differentiation  with  respect  to  n 

Denote  differentiation  with  respect  to  £ 
Denotes  a  vector  quantity 


Superscript 

*  Denotes  a  quantity  in  the  transformed  plane 

-  Denotes  a  turbulence  time  averaged  quantity 

1  Denotes  a  turbulent  fluctuating  quantity 

Denotes  a  temporary  nondimensional  variable  definition 

s  Iteration  number  in  the  successive-over-relaxation 

iteration  procedure 

n  Pertaining  to  n  contours 

£  Pertaining  to  £  contours 


Prefix 

d 

D 

A 

l 

V 


Derivative  operator 

Substantial  derivative  operator  (D/Dt  *  a -/a t  +  v/V-) 

Increment 

Summation 

Del -operator  defined  as  i  —  +  j— 
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SECTION  I 
INTRODUCTION 

Much  effort  has  been  expended  by  the  aeronautical  community  in 
determining  the  aerodynamic  characteristics  of  airfoils.  Linear  methods 
are  extensively  used  in  design  work  for  studying  configurations  at  small 
angles  of  attack  with  negligible  flow  separation.  Experimental  wind 
tunnel  investigations  are  used  to  determine  the  characteristics  near  stall 
where  separation  phenomena  become  important.  Recent  developments  in 
numerical  techniques  have  stimulated  research  on  another  approach,  namely 
the  numerical  solution  of  the  Navier-Stokes  equations.  These  equations 
model  the  viscous  effects  which  contribute  to  airfoil  stall.  For  this 
reason,  numerical  Navier-Stokes  methods  offer  the  possibility  of  determin¬ 
ing  the  aerodynamic  characteristics  for  airfoils  experiencing  stall. 
Numerical  methods  can  also  complement  experimental  methods  by  efficiently 
extending  the  range  of  parameters  under  investigation.  Furthermore, 
numerical  methods  eliminate  model  support  interference  and  wall  interference 
effects  found  in  wind  tunnel  testing. 

The  purpose  of  this  investigation  is  to  develop  a  numerical  Navier- 
Stokes  method  that  will  accurately  determine  the  aerodynamic  characteristics 
of  incompressible  turbulent  viscous  flow  over  two-dimensional  airfoils  near 
stall . 

The  development  of  a  numerical  method  for  turbulent  flow  requires 
a  suitable  technique  for  distributing  points  throughout  the  flow  field 
and  a  model  which  describes  the  behavior  of  turbulence  within  specific 
regions  of  the  flow  field.  A  survey  of  numerical  grid  generating  techniques 
and  available  turbulence  models  is  presented.  The  quantity  of  literature 
concerned  with  the  Navier-Stokes  equations  is  extensive.  Therefore  a 
summary  of  the  literature,  which  describes  formulations  of  the  Navier- 
Stokes  equations  and  their  numerical  solving  techniques  applied  to  flows 
over  airfoils,  is  given.  The  research  objectives  for  this  work  are 
then  discussed,  and  a  summary  of  the  remaining  sections  is  presented. 
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Body-fitted  curvilinear  coordinate  systems  greatly  enhance  the 
application  of  numerical  methods  to  practical  boundary  value  problems 
involving  partial  differential  equations.  The  representation  of  a 
boundary  surface  as  a  coordinate  line  reduces  the  difficulties  associated 
with  numerically  specifying  boundary  conditions  by  interpolation  in  finite 
difference  methods.  In  the  physical  (x,y)  plane,  values  for  one  compu¬ 
tational  coordinate  ?  are  specified  at  selected  points  on  both  the  body 
surface  and  the  outer  boundary.  Constant  values  for  the  other  computational 
coordinate  n  are  specified  on  both  the  body  surface  and  the  outer  flow 
boundary.  The  transformed  computational  (c,n)  plane  then  becomes  a 
rectangular  region  with  an  orthogonal  grid.  Winslow  (Reference  1)  and 
Chu  (Reference  2)  introduced  the  concept  for  two-dimensional  regions 
interior  to  a  closed  boundary.  Their  transformed  coordinates  are  solutions 
of  Laplace's  equation  in  the  physical  plane  and  define  a  triangular  mesh 
system  in  the  physical  plane.  Amsden  and  Hirt  (Reference  3)  took  the 
physical  coordinates  to  be  solutions  of  a  modified  Laplace's  equation 
in  the  transformed  plane. 

Thompson,  Thames,  and  Mastin  (References  4,5,6)  generalized  the 
method  for  the  automatic  generation  of  body-fitted  coordinates  for  any 
two-dimensional,  multi -connected  region.  They  also  introduced  the  use 
of  forcing  functions  in  a  Poisson  equation  which  provides  mesh  control 
in  regions  with  large  gradients.  Hodge  (Reference  7)  developed  an 
automated  grid  line  attraction  method  based  on  boundary-layer  theory 
which  determines  the  coefficients  in  his  forcing  function.  Hodge 
assumed  a  Blasius  boundary-layer  profile  and  distributed  his  grid 
points  at  approximately  equal  velocity  increments  in  the  boundary-layer. 
Steger  and  Sorenson  (Reference  8)  introduced  auxiliary  conditions  for 
the  forcing  functions  which  provide  angle  and  distance  control  at  the 
inner  boundary  surface.  The  angle  with  which  a  £  line  intersects  the 
body  surface  is  specified  by  a  function  0(5),  and  the  rate  of  change 
of  arclength  with  n  on  a  £  line  at  the  body  surface  is  prescribed  by 
s^U).  Sorenson  (Reference  9)  later  imposed  similar  conditions  on 

an  outer  computational  boundary.  These  auxiliary  conditions  are  used  to 
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solve  for  coefficients  in  the  chosen  exponential  forcing  functions.  The 
geometric  conditions  hold  exactly  only  in  the  limit  as  An  approaches  zero. 
Sorenson  reported  that  numerical  instabilities  occur  for  large  changes 
in  the  coefficients  during  successive  iterations  and  for  boundaries  with 
sharp  corners.  He  implemented  a  limit  function  which  damped  the  change  of 
the  value  for  each  coefficient  over  one  iteration;  and  at  sharp  corners  ne 
computed  average  values  of  each  coefficient  from  data  at  the  neighboring 
boundary  points.  Mastin  and  Thompson  (Reference  10)  have  also  extended 
the  elliptic  body-fitted  coordinate  generation  technique  to  three 
dimensions  for  simple  geometries. 

Some  other  specific  grid  generation  techniques  which  use  properties 
of  elliptic  differential  equations  have  also  been  reported.  Meyder 
(Reference  11)  constructed  an  orthogonal  curvilinear  coordinate  system 
by  using  electric  field  theory.  He  solved  the  potential  equation  twice 
with  different  sets  of  mixed  Dirichlet  and  Neumann  boundary  conditions 
for  the  electric  potential  and  electric  force  lines,  respectively. 

These  solutions,  however,  were  obtained  in  the  physical  plane  using 
interpolation  and  became  the  curvilinear  coordinate  lines  in  the 
physical  plane.  The  coordinate  metrics  were  then  used  to  formulate  a 
finite  difference  equation  which  was  also  solved  in  the  physical  plane. 
Conformal  mapping  techniques  which  employ  the  Theodorsen-Garrick 
(Reference  12)  transformation  have  been  examined  by  Ives  (Reference  13). 

He  introduced  the  use  of  fast  Fourier  transform  methods  and  developed  a 
new  class  of  transformations  which  map  the  flow  field  of  a  two-element 
airfoil  onto  the  region  between  two  concentric  circles.  Conformal 
techniques  are  not  capable  of  being  extended  to  three-dimensional 
geometries  (Reference  4).  Neither  of  these  approaches  offers  a 
convenient  means  for  mesh  control  in  regions  of  large  flow  gradients. 

Recently,  a  geometric  grid  generation  technique  has  been  introduced 
by  Glbeling,  Shamroth,  and  Eiseman  (Reference  14).  The  technique 
parameterizes  (t)  the  body  and  outer  boundary  surfaces  and  uses  a 
stretching  function  R(r)  for  mesh  attraction  along  constructed  lines 
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perpendicular  to  the  body  surface.  Unit  increments  for  ordered  pairs  (t,r) 
generate  the  corresponding  computational  plane.  Further  refinements  which 
provide  angle  and  arclength  variation  control  at  a  body  surface  have 
subsequently  been  developed  by  Eiseman  (References  15,16,17). 

The  search  for  an  accurate  and  universal  turbulence  model  has 

paralleled  the  development  of  body-fitted  grid  generation  techniques. 

In  principle,  the  Navier-Stokes  equations  completely  describe  the 

turbulent  fluctuating  fluid  motion.  The  required  mesh  resolution, 

necessary  to  resolve  the  turbulent  eddies  with  varying  length  scales, 

when  translated  into  computer  resources  presently  make  this  approach 

unfeasible.  Many  quantities  of  engineering  interest  in  turbulent  flows 

involve  a  mean  value  taken  over  a  time  interval.  The  time  interval  is 

sufficiently  long  to  include  many  fluctuations  while  small  compared  with 

the  characteristic  time  of  the  mean  flow.  The  Navier-Stokes  equations 

can  then  be  reformulated  using  these  mean  flow  variables.  This  Reynolds 

averaging  procedure  introduces  the  mean  of  terms  involving  fluctuating 

quantities.  The  Reynolds  stress  components  -puTuT  are  the  most  common 

*  J 

terms  of  this  type.  To  solve  the  averaged  form  of  the  Navier-Stokes 
equations,  "turbulence  closure"  must  be  achieved  by  suitably  modelling 
these  additional  terms.  This  approach  to  turbulence  has  led  to  models 
which  introduce  auxiliary  relationships  ranging  from  algebraic  equations 
to  several  partial  differential  equations.  These  models  are  commonly 
categorized  and  are  now  summarized  according  to  the  number  of  additional 
partial  differential  equations  which  comprise  the  model. 

The  algebraic  or  zero  equation  models  have  their  origins  in 
Boussinesq's  (Reference  18)  eddy  viscosity  hypothesis  and  Prandtl’s 
(Reference  19)  mixing  length  model.  The  local  turbulent  stresses  are 
assumed  proportional  to  the  local  mean  flow  strain  rates  with  the 
proportionality  constant  defined  as  an  equivalent  or  eddy  viscosity. 

The  eddy  viscosity  models  of  Cebeci-Smith  (Reference  20),  Mellor-Herring 
(Reference  21),  and  Patankar-Spalding  (Reference  22)  represent  this 
approach.  The  boundary-layer  in  each  method  is  divided  into  an  inner 
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near  wall  layer  and  an  outer  wake  layer  with  separate  expressions  for  the 
eddy  viscosity  coefficient.  In  the  Cebeci-Smith  and  Patanker-Spalding 
models  the  inner  mixing  length  varies  as  a  linear  function  of  the  normal 
distance  from  the  wall  modified  with  a  Van  Driest  (Reference  23)  laminar 
sublayer  correction.  The  Cebeci  and  Mellor  models  both  base  the  outer 
length  scale  on  the  displacement  thickness  while  the  Patanker-Spalding 
model  uses  the  boundary- layer  thickness  directly.  These  models  have  been 
successfully  extended  and  applied  (Reference  24)  to  a  wide  variety  of 
boundary-layer  flow  geometries  involving  compressibility,  heat  and  mass 
transfer,  and  curvature  effects.  In  addition,  the  mean  flow  models  are 
computationally  efficient.  Launder  and  Spalding  (Reference  25)  point 
out  that  these  models  predict  a  vanishing  eddy  viscosity  where  the 
velocity  gradient  is  zero  and  have  not  been  successful  for  large 
separated  recirculating  flows. 


The  one  equation  turbulence  model  introduced  by  Prandtl  (Reference 
26)  is  an  extension  of  the  algebraic  technique.  In  this  approach  the 
solution  of  the  partial  differential  turbulent  kinetic  energy  equation 
usually  provides  a  local  velocity  scale  given  by  q2  =  u^ul  where 

summation  is  implied.  The  length  scale  L  is  prescribed  by  an  algebraic 
expression  as  before.  Terms  involving  fluctuating  quantities  other  than 
q  must  still  be  modelled.  Glushko  (Reference  27),  Mellor  and  Herring 
(Reference  28),  and  Wolfshtein  (Reference  29)  model  the  gradient  diffusion 
term  -((u1^  p/p  )  +  u'^  (u1^  U'.  )/2)  with  a  gradient  of  q2  given  as 
2 

NQ  eM  9k  ^  where  Ng  is  a  specified  constant  or  function  and  is 

the  eddy  viscosity.  The  Reynolds  stress  is  related  to  the  mean  flow  as 
in  the  algebraic  approach  except  that  the  eddy  viscosity  is  assumed 
proportional  to  qL.  Bradshaw,  et  al  (Reference  30)  model  the  gradient 
diffusion  term  with  an  expression  G  q2  where  G  is  an  empirical 


constant  or  function  and  is  a  velocity  characteristic  of  large 
eddy  motions.  They  also-  write  the  Reynolds  stress  directly  as  a 


function  of  q  in  the  form  x 


tij 


A.  .  q 
ij 


In  each  model  the  isotropic 


dissipation  is  modelled  by  a  form  q3/L.  Nee  and  Kovasznay  (Reference  31) 
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use  a  rate  equation  for  the  total  viscosity  N  =  u  +  in  place  of  the 

kinetic  energy  equation.  They  formulated  expressions  for  the  generation 
and  dissipation  terms  involving  various  constants  and  the  length  scale. 
Launder  and  Spalding  (Reference  25)  observe  that  the  one  differential 
equation  models  require  a  moderate  increase  in  computer  resources  but  do 
not  in  general  provide  more  accurate  results  than  the  results  obtained 
from  algebraic  methods. 

Kolmogorov  (Reference  32)  introduced  the  somewhat  more  complex  type 
of  turbulence  model  which  uses  two  partial  differential  equations.  A 
form  of  the  turbulent  kinetic  energy  equation  provides  a  local  velocity 
scale.  The  local  length  scale  is  obtained  from  a  second  equation.  Ng 
and  Spalding  (Reference  33)  formalized  this  approach  by  introducing  the 
energy-length  equation  derived  by  Rotta  (Reference  34).  They  also  used 
the  Glushko  closure  model  in  the  turbulent  energy  equation.  This  class 
of  two  equation  models  is  named  the  k  -  kL  turbulence  model  where  the 
eddy  viscosity  is  proportional  to  k^  L.  Saffman  (Reference  35)  has  used 
a  transport  equation  for  the  mean  vorticity  w  together  with  the  turbulent 
kinetic  energy.  The  eddy  viscosity  becomes  proportional  to  k/w  in  this 
k  -  o)  class  of  closure  models  for  turbulence.  The  local  length  scale 
is  assumed  to  be  proportional  to  k2/w.  Wilcox  and  Rubesin  (Reference 
36)  have  modified  this  approach  for  compressible  flows  and  generalized 
the  constitutive  equation.  Jones  and  Launder  (Reference  37)  use  a 
transport  equation  for  the  rate  of  dissipation  of  turbulent  energy  e 
along  with  the  turbulent  kinetic  energy  equation.  Glushko  type  closure 
models  are  assumed  for  the  terms  involving  fluctuating  quantities  in 
both  equations.  This  k  -  e  class  of  two  equation  models  has  a  local 

3/2 

length  scale  proportional  to  k  /e  with  the  eddy  viscosity  proportional 

2 

to  k  /e.  The  two  equation  models  have  been  applied  to  various  boundary- 
layer  and  free  shear  layer  flows  with  a  variety  of  constants  and  closure 
models.  The  two  equation  models  require  significant  increases  in 
computational  resources  and  have  not  led  to  a  model  of  universal 
applicability  (References  24,25). 
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The  search  for  a  more  general  turbulence  model  has  led  to  the  use  of 
the  transport  equations  for  the  Reynolds  stresses.  In  this  approach  the 
eddy  viscosity  concept  is  discarded.  However,  closure  models  are  still 
required  for  the  terms  containing  fluctuating  quantities  other  than  the 
Reynolds  stresses.  Donaldson  (Reference  38)  introduces  closure  models 
which  express  the  fluctuating  terms  as  functions  of  the  Reynolds  stress 
and  chosen  length  scales.  Hanjalic  and  Launder  (Reference  39)  introduce 
closure  models  which  use  the  Reynolds  stress  and  retain  the  equations  for 
turbulent  kinetic  energy  k  and  turbulent  dissipation  e.  Only  cases  with 
the  simplified  two-dimensional  boundary- layer  approximations  have  been 
investigated.  The  extensive  computational  resources  and  initial  state 
of  development  of  these  models  preclude  this  type  of  approach  from  practical 
consideration  as  a  turbulence  model  for  a  complex  flow  field  computation. 

Theoretical  turbulence  models  of  further  complexity  have  appeared. 
Kolovandin  and  Votutin  (Reference  40)  introduced  a  statistical  theory 
where  additional  equations  are  obtained  for  other  correlations  of  the 
fluctuating  quantities.  Ferziger  (Reference  41)  took  a  meteorological 
viewpoint  of  turbulent  flows  by  numerically  simulating  large  scale  eddies 
while  modeling  the  small  scale  structures  with  an  eddy  viscosity  technique. 
This  approach  is  a  first  step  toward  numerically  solving  a  turbulent  flow 
field  with  the  instantaneous  equations  of  motion. 

The  review  of  available  turbulence  models  provides  the  basis  for 
selecting  a  suitable  approach  for  use  in  the  numerical  solution.  The 
algebraic  and  the  two  equation  eddy  viscosity  models  appear  to  be 
realistic  choices.  As  previously  reported  (Reference  25),  the  one 
equation  techniques  yield  unimproved  results  compared  with  algebraic 
methods  and  at  additional  cost  in  computer  resources.  The  two  equation 
methods  require  the  solution  of  additional  partial  differential  equations. 

In  these  methods,  terms  involving  fluctuating  quantities,  except  the 
Reynolds  stresses,  must  still  be  modelled  using  additional  coefficients. 

For  these  reasons,  the  simpler  algebraic  eddy  viscosity  technique,  which 
requires  considerably  less  computer  resources,  is  employed.  If  the 
physical  phenomena  associated  with  separated  adverse  pressure  gradient 
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flows  can  be  included  in  the  turbulence  model,  then  the  accurate 
calculation  of  the  aerodynamic  characteristics  may  be  accomplished 
with  an  algebraic  technique. 

A  survey  of  the  previous  research  involving  numerical  solutions  of 
the  Navier-Stokes  equations  for  flow  over  airfoils  establishes  the  pre¬ 
diction  level  of  computational  methods  and  also  reviews  the  numerical 
algorithms.  The  grid  generation  technique  and  turbulence  model  employed, 
when  applicable,  are  also  included. 

Early  numerical  solutions  of  the  Navier-Stokes  equations  for  flow  over 
airfoils  used  the  vorticity-stream  function  formulation  with  an  automated 
grid  generation  technique  (Reference  4).  Walker  (Reference  42)  applied 
the  method  to  the  laminar  flow  over  a  flat  plate  and  compared  the  numerical 
solution  with  Blasius*  (Reference  43)  solution.  Thames  (Reference  44) 
used  body-fitted  coordinates  with  the  vorticity-stream  function  approach 
and  solved  the  Navier-Stokes  equations  for  various  bodies  in  doubly 
connected  regions.  He  obtained  solutions  for  the  flow  over  airfoils  at 

4 

chord  Reynolds  numbers  less  than  10  .  Problems  mainly  attributed  to  wall 
vorticity  developed  for  solutions  of  airfoils  at  angle  of  attack. 

Mehta  and  Lavan  (Reference  45)  also  used  a  vorticity-stream  function 
method  in  studying  the  laminar  starting  vortex  and  separation  bubbles  for 
impulsively  started  incompressible  laminar  flow  over  a  Joukowski  airfoil 
at  Reynolds  numbers  less  than  10^.  They  used  three  point  backward  time 
differences  and  centered  spatial  differences.  The  computational  grid 
was  obtained  through  a  conformal  transformation  followed  by  a  radial 
stretching  transformation. 

Reddy  and  Thompson  (Reference  46)  applied  an  integro-differential , 
vorticity-velocity  field  method  for  the  solution  of  incompressible  flow 
in  doubly  connected  regions.  Backward  time,  centered  spatial  (BTCS) 
differences  were  applied  to  the  Navier-Stokes  equations.  The  difference 
equations  were  solved  using  successive-over-relaxation  (SOR)  iteration. 
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They  also  employed  the  coordinate  mesh  attraction  technique  with  a  time 
dependent  expanding  mesh  system.  Symmetric  airfoils  at  zero  angle  of 
attack  with  a  Reynolds  number  less  than  10^  were  considered.  At  the 
higher  Reynolds  numbers,  the  calculation  of  surface  vorticity  required 
a  large  number  of  grid  points  on  the  surface  which  greatly  increased  the 
required  computer  time.  A  steady  state  solution  was  not  obtained. 

The  vorticity-velocity  field  formulation  has  been  applied  by  Sankar 

and  Wu  (Reference  47)  to  the  case  of  incompressible  laminar  flow  about  an 

oscillating  airfoil.  They  used  a  12%  thick  Joukowski  airfoil  at  a  Reynolds 

number  of  1000.  Triangular  finite  elements  were  constructed  near  the 

airfoil  surface  and  a  rectangular  mesh  was  used  away  from  the  airfoil . 

A  conformal  transformation  was  used  to  map  the  airfoil  in  the  physical 

plane  onto  a  unit  circle  in  the  computational  plane.  Sankar  and  Tassa 

(Reference  48)  investigated  this  same  problem  with  a  compressible  flow 

formulation  of  the  Navier-Stokes  equations.  They  used  a  conformal 

transformation  followed  by  an  algebraic  radial  stretching  transformation. 

The  alternating-direction-implicit  (ADI)  finite  difference  method  of 

Briley  and  McDonald  (Reference  49)  was  used  to  obtain  solutions  for 

a 

Reynolds  numbers  less  than  10  and  a  Mach  number  of  0.2.  In  a  separate 
research  effort,  Sugavanam  and  Mu  (Reference  50)  attempted  to  use  a  two 
equation  k-e  turbulence  model  with  a  vorticity-velocity  formulation. 

A  conformal  transformation  for  a  12  percent  thick  Joukowski  airfoil 
was  used.  They  experienced  difficulties  in  obtaining  a  converged 
solution  even  after  eight  hours  of  CYBER  74  CPU  time  were  expended. 
Variations  for  both  lift  and  drag  of  the  order  of  50  percent  occurred. 

The  use  of  the  primitive  variables  of  velocity  and  pressure  for 
incompressible  flow  was  introduced  by  Harlow  and  Welch  (Reference  51) 
in  the  explicit  forward  time,  centered  spatial  (FTCS)  Mark-and-Cell  (MAC) 
method.  They  included  a  Poisson  equation  for  pressure  which  is  obtained 
by  taking  the  divergence  of  the  momentum  equation.  They  also  found  that 
the  velocity  divergence  terms  in  the  momentum  equations  were  required 
for  the  pressure  field  calculation.  Hirt  and  Harlow  (Reference  52) 
further  developed  the  method.  Hodge  (Reference  53)  considered  the  case 
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of  laminar  incompressible  viscous  flow  for  an  airfoil  at  angle  of  attack. 
He  used  a  form  of  the  Thompson,  et  al  (Reference  4)  body-fitted  grid 
transformation  and  applied  an  implicit  BTCS  differencing  method  to  the 
Navier-Stokes  equations.  The  system  of  difference  equations  was  solved 
with  SOR  iteration.  Various  methods  of  calculating  the  pressure  field 
were  investigated.  Hodge  concluded  that  the  momentum  equations  should 
retain  the  velocity  divergence  terms  and  that  a  Poisson  pressure  equation 
should  be  used  to  satisfy  continuity.  Ghia,  Hankey,  and  Hodge  (Reference 
54)  applied  a  Poisson  pressure  equation  and  the  primitive  variables  form 
of  the  Navier-Stokes  equations  to  study  incompressible  driven  flow  in  a 
square  cavity  for  Reynolds  numbers  under  1000.  They  used  an  alternating- 
direction-implicit  (ADI)  finite  difference  technique  for  the  momentum 
equations  and  SOR  iteration  for  the  pressure  equation.  A  Neumann  boundary 
condition  derived  from  the  normal  component  of  the  momentum  equations  was 
employed  to  compute  the  wall  pressure. 

The  Poisson  pressure  equation  has  been  further  examined  by  Chien 
(Reference  55)  for  internal  flows.  He  observed  that  the  calculation  of 
pressure  by  direct  integration  of  the  momentum  equations  can  be  inaccurate 
when  large  velocity  variations  are  present.  He  found  that  forms  of  the 
pressure  equation  were  more  suitable  for  computing  the  pressure  field 
near  boundaries.  Recently,  Hodge,  et  al  (Reference  7)  applied  an  implicit 
backward  time  finite  difference  method  with  both  upwind  and  centered 
spatial  differences  to  the  Navier-Stokes  equations  in  primitive  variables 
form.  A  Poisson  pressure  equation  was  used  to  obtain  a  solution  for 
laminar  flow  over  airfoils  at  angle  of  attack  with  a  Reynolds  number  of 

5 

order  10  .  The  solution  contains  a  large  oscillating  separated  region 
which  gives  approximate  trends  in  lift  but  very  poor  agreement  with 
available  drag  data  (Reference  56). 

Early  numerical  Navier-Stokes  investigations  of  turbulent  flow  over 
airfoils  considered  two-dimensional  transonic  viscous  flow  at  small  angles 
of  attack.  Delwert  (Reference  57)  applied  MacCormack's  explicit  method 
(Reference  58)  with  a  50x38  rectangular  based  exponentially  stretched 
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mesh.  He  used  20  uniformly  spaced  points  to  define  the  upper  airfoil 
surface  and  imposed  the  Neumann  boundary-layer  condition  3p/3n  =  0  on  the 
surface.  Steger  (Reference  59)  used  the  implicit  Beam-Warming  method 
(Reference  60)  with  a  71x33  grid  obtained  from  a  modified  Thompson 
transformation  (Reference  4).  He  calculated  the  solution  on  the 
branch  cut  with  a  linear  extrapolation  procedure  and  evaluated  surface 
pressures  using  the  momentum  equations  for  the  direction  normal  to  the 
surface.  Walitt,  et  al  (Reference  61)  applied  the  first  order  method  of 
Trulio  (Reference  62)  coupled  with  a  boundary-layer  technique  used  near 
the  airfoil  leading  edge.  They  used  a  130x68  mesh  system  but  did  not 
capture  the  full  suction  pressure  on  the  upper  surface  near  the  airfoil 
nose.  Each  investigation  used  a  Cebeci -Smith  (Reference  63)  type 
algebraic  eddy  viscosity  model  and  the  Reynolds  averaged  compressible 
Navier-Stokes  equations. 

Recently,  Shamroth  and  Gibeling  (Reference  64)  used  the  Briley- 
McDonald  (Reference  49)  implicit  finite  difference  formulation  with  a 
constructive  type  81x30  grid  system  to  compute  turbulent  flow  over  an 
airfoil  at  an  angle  of  attack  of  six  degrees  and  Reynolds  number  of  106. 
They  first  tried  a  two  equation  k-e  turbulence  model  but  experienced 
convergence  problems  near  the  leading  edge  and  far  field  flow  regions. 
They  reported  obtaining  large  or  negative  values  for  the  turbulent 
viscosity  in  essentially  laminar  regions.  The  turbulent  kinetic  energy 
equation  with  an  algebraical ly  prescribed  length  scale  was  then  used. 
Major  discrepancies  occurred  for  the  mean  pressure  distribution  solution 
in  the  suction  peak  and  trailing  edge  regions  of  the  airfoil  surface. 

They  also  noted  that  the  constructive  grid  technique  caused  a  crossing 
of  grid  lines  of  the  same  family  for  highly  cambered  airfoils. 

The  examination  of  existing  numerical  solutions  of  the  Navier-Stokes 
equations  for  two-dimensional  fow  over  airfoils  reveals  Reynolds  number 
and  angle  of  attack  restrictions.  As  a  result,  numerical  Navier-Stokes 
methods  for  the  accurate  computation  of  the  flow  field  and  resulting 
aerodynamic  characteristics  for  an  airfoil  near  stall  are  unavailable. 
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The  development  of  such  a  numerical  method  for  solving  incompressible 
two-dimensional  turbulent  flow  over  airfoil  sections  near  stall  does 
constitute  a  significant  contribution  in  computational  aerodynamics. 

The  fulfillment  of  this  research  goal  requires  the  formulation  of  an 
adequate  turbulence  model  for  use  in  both  leading  and  trailing  edge 
separated  regions.  The  far  field  boundary  conditions  for  incompressible 
flows  must  also  be  examined. 

The  selection  of  a  suitable  numerical  approach  for  incompressible 
turbulent  flows  is  based  on  several  considerations.  An  implicit  technique 
is  preferred  because  of  the  required  small  grid  spacing  near  the  body 
surface  necessary  to  resolve  viscous  stresses.  Stability  criteria  of 
explicit  methods  would  impose  a  small  time  step  restriction  as  a  result 
of  the  grid  size.  The  numerical  method  should  be  capable  of  predicting 
laminar  flow  for  Reynolds  numbers  approaching  turbulent  flow  conditions 
since  regions  of  laminar  flow  may  occur.  The  use  of  primitive  variables 
is  readily  extended  to  three  dimensions  and  provides  for  the  direct  calcu¬ 
lation  of  the  pressure  field.  These  criteria  are  satisfied  by  an  implicit 
finite  difference  procedure  developed  by  Hodge  (References  53,7)  which 
uses  primitive  variables  and  successive-over-relaxation  iteration.  The 
implicit  finite  difference  method  is  employed  in  this  investigation. 

The  remaining  sections  of  this  research  effort  discuss  the  techniques 
employed  in  achieving  the  defined  research  goal.  In  Section  II  a  form  of 
the  Thompson  numerically  generated  body-fitted  coordinate  system  is 
described.  The  methods  of  attracting  grid  lines  to  the  body  surface 
and  distributing  surface  points  are  also  discussed.  The  governing 
equations  for  incompressible  two-dimensional  turbulent  flow  in  primitive 
variables  are  presented  in  Section  III.  Boundary  and  initial  conditions 
are  discussed.  An  algebraic  eddy  viscosity  turbulence  model  formulated 
for  adverse  pressure  gradient  separated  flows  is  described.  An  algebraic 
type  of  turbulence  model  has  been  chosen  because  of  the  previous 
computational  successes  for  other  turbulent  flow  problems  and  the  definite 
computational  efficiency  obtained  by  this  approach.  The  method  for 
evaluating  force  coefficients  exerted  on  the  body  by  the  flow  field  is 
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then  discussed.  The  numerical  techniques  used  in  obtaining  a  solution 
are  presented  in  Section  IV.  First,  the  finite  difference  method  for 
numerically  generating  the  grid  is  given.  The  Navier-Stokes  equations 
are  then  written  in  finite  difference  form.  The  numerical  method  used 
to  implement  the  boundary  and  initial  conditions  is  described.  The 
computational  procedure  for  introducing  turbulence  is  discussed  next. 

The  finite  difference  and  numerical  integration  techniques  for  calculating 
the  force  coefficients  are  then  presented.  The  formal  discussion  of  the 
numerical  solutions  and  comparisons  with  experimental  data  are  given  in 
Section  V.  The  conclusions  derived  from  this  research  work  and 
recommendations  for  future  work  are  set  forth  in  Section  VI. 
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SECTION  II 

BODY-FITTED  CURVILINEAR  COORDINATE  SYSTEM 

In  this  section  the  technique  is  described  which  numerically  generates 
body-fitted  curvilinear  coordinate  systems  as  previously  developed  by 
Thompson,  Thames,  and  Mastin  (References  4,5).  Coordinate  lines  are 
attracted  to  the  body  surface  with  a  method  introduced  by  Hodge 
(Reference  7).  The  distribution  of  grid  points  which  define  the 
body  surface  is  discussed. 

1 .  TRANSFORMATION 

The  doubly-connected,  two-dimensional  region  R  bounded  by  closed 
contours  Cb  and  C^.  of  arbitrary  shape  is  to  be  transformed  into  a 

rectangular  domain  R*  as  shown  in 
from  the  physical  (x,  y)  plane  to 
plane  is  given  by 

£  = 

The  inverse  transformation,  if  it 

x  = 

The  related  Jacobian  is  given  as 

J  = 

where  the  subscripts  denote  differentiation. 

The  existence  of  the  inverse  transformation  has  been  studied  by 
Thames  (Reference  44).  He  showed  that  if  the  functions  of  Equation  1 
are  continuously  differentiable  at  a  point  and  the  Jacobian  is  nonsingular, 
then  the  general  transformation  has  an  Inverse  in  a  disk  about  the  point. 


Figure  1.  The  general  transformation 
the  transformed  or  computational  (£,n) 


£  (x,y)  n  =  n  (x,y) 


(1) 


exists,  can  be  expressed  as 


x{£,  n)  y  =  y(£,  n) 


(2) 


X£  yn  '  Xn  y£ 


(3) 
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In  addition,  if  the  £  and  n  transformation  functions  are  harmonic,  then 
a  weak  maximum  principle  is  obeyed.  Then,  no  extrema  will  exist  in  region 
R  and  consequently  the  Jacobian  J  will  not  be  zero  there. 

Partial  derivatives  are  transformed  using  the  chain  rule  and  inverse 
transformation  relations.  For  a  sufficiently  differentiable  function  f 
of  x  and  y,  the  derivatives  transform  as  follows: 


fx  ■  <*„  V/J 

(4) 

f  =  (-x  f  +  X.  f  I/O 

(5) 

y  n  t,  i,  n 

Higher  order  derivatives  are  similarly  derived. 

The  system  of  elliptic  coordinate  generating  equations 

is  chosen 

to  be 

Exx  *  Eyy  ■  p(!-  "> 

(6a) 

"xx  *  V  '  0<E'  n) 

(6b) 

with  the  boundary  conditions  on  curve 

£  =  £b(*>  y)  n  =  nb  (7) 


and  on  curve  C^. 


£  —  £^(x,  y)  n  n^ 


(8) 


where  £b(x,  y)  and  £^(x,  y)  are  prescribed  functions  and  and 
are  prescribed  constants.  The  other  two  sides  C*  and  c£  in  the 
rectangular  region  are  transformed  from  the  branch  cut  -  Cg. 

The  functions  x(c,  n)  and  y(s,  n)  and  all  derivatives  are  continuous 
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across  this  cut.  These  boundary  conditions  ensure  that  the  body  and  the 
far  field  boundaries  are  defined  by  constant  n  lines.  The  generalized 
functions  P  and  Q  are  utilized  in  attracting  the  coordinate  lines  in 
various  regions. 

The  rectangular  grid  system  in  the  transformed  plane  is  used  for 
the  computations.  The  mesh  generating  system  (Equations  6a  and  6b)  is 
transformed  to  the  computational  plane  and  becomes 


where 


55 


“  •»  1  Y  A 

Cn  nn 


'5n  nn 


[x?P(s,  n) 

+  x  Q(e,  n)] 

(9a) 

[y^P(5.  n) 

+  ynQ(C,  n)] 

(9b) 

(10) 

(ID 

:  +  ycy 

n  5  n 

(12) 

The  boundary  conditions  (Equations  7  and  8)  similarly  become  on  curve  C. 


x  =  fb(5,  nb)  y  =  gK(5,  nK) 


bv<”  'V 


(13) 


and  on  curve  C. 


x  ~  f^(5f  n.^)  ¥  ~  9f ^ >  n.^) 


(14) 


The  functions  fb>  f^.,  gb>  and  g^  are  determined  by  the  known  contours  Cb 

f 


and  Cf  and  the  functions  Cb  and  e; 
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2.  COORDINATE  LINE  ATTRACTION 

The  resolution  of  viscous  terms  near  the  body  surface  requires  a 
fine  mesh  spacing  in  this  region.  The  forcing  functions  P  and  Q  are  used 
to  attract  £  and  n  lines,  respectively.  In  this  work  the  £  line  spacing 
is  accomplished  by  specifying  the  boundary  condition  function  ^(x,  y). 

The  n  line  spacing  is  adjusted  through  the  use  of  the  Q  function  proposed 
by  Thompson  (Reference  4)  and  modified  by  Hodge  (Reference  7). 

K(£) 

Q(£,  n)  =  -  £  Aj£)  exp[-D(?)|n  -  nj]  (15) 
k=l  K  K 


where  Ak  are  the  amplification  factors,  D  is  the  damping  factor,  K  is  the 

number  of  terms,  and  nk  =  (k  -  l)An.  The  amplification  and  damping  factors 

are  determined  by  requiring  that  the  tangential  velocity  u  be  a  linear 
function  of  n.  In  this  case,  the  second  and  higher  order  derivatives 
with  respect  to  n  become  zero  thereby  minimizing  the  truncation  error 
of  u  in  the  transformed  plane. 


Hodge  (Reference  7)  has  developed  a  technique  which  approximately 
satisfies  the  above  criterion.  The  required  point  distribution  on  a  £ 
line  is  determined  using  Blasius'  fiat  plate  boundary-layer  solution. 

The  similarity  series  solution  (Reference  43)  is  given  by 

00  C  a  n+1 

f  (a )  =  l  {-hf  J1—! -  a3n+2  (16) 

(3n  +  2)  ! 


where  the  similarity  variable  a  =  y[u  /vx] 2.  u.  is  the  edge  boundary- 

S  6 

layer  velocity  and  v  is  the  kinematic  viscosity.  Differentiating  Equation 
16  with  respect  to  a  gives  an  expression  for  the  tangential  velocity  u  in 
the  boundary-layer  , 


f '  (a )  =  = 

ue 


(-h)n 


C  a. 
n  1 


3n+1 


n=0 


(3n  +  1)! 


(17) 
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The  first  four  coefficients  of  the  Blasius  series  are  Cq  =  1 ,  =1, 

C2  =  11,  and  =  375.  Also,  from  Equation  17  the  factor  is  related 
to  f  by  =  f"(0)  which  is  0.33206.  Equation  17  is  conveniently  written 

as  f‘(o)  -  jp  =  0  (18) 

e 


Then,  for  a  given  £  line  (which  determines  x)  and  a  selected  number  of 
points  in  the  boundary-layer,  values  for  a  which  give  equal  tangential 
velocity  increments  are  calculated  from  Equation  18.  The  corresponding 
y  distances  are  computed  by  using  the  definition  of  o  which  yields 


y 


(19) 


These  computed  y  distances  become  the  desired  points  near  the  body  for  the 
n  line  intersections  with  the  specified  £  line. 


The  second  part  of  the  boundary-layer  grid  line  attraction  technique 
involves  calculating  values  of  Ak  and  D  which  when  used  in  Q  give  the 

desired  n  line  spacing.  Near  the  body  surface,  the  n  transformation 
equation  (Equation  6b)  can  be  written 


(O  = 

yyb 


Q(£,  n) 


(20) 


where  variations  with  respect  to  x  are  neglected.  If  the  expression  for 
the  forcing  function  Q  in  Equation  15  is  substituted  into  Equation  20, 
then  Equation  20  can  be  integrated  with  respect  to  n  to  obtain 


K 

s  2 ( A^/D)  {exp[-D  H(n,  nk)]  +  H(sgn(nk  -  n),  0.) 


[1  -  exp(-D)n  -  nkl )]) 


(21) 


where  the  difference  function  H  and  sgn  function  are  defined  as  follows 
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0 

H(n,  nk)  = 

if 

n  <  nk 

n  -  nk 

if 

n  >  nk 

1 

if 

n.  -  n  >  0 

sgn(nk  -  n)  = 

K  — 

-1 

if 

35 

A 

O 

Now  evaluate  the  expression  for  (dn/dy;  (Equation  21)  at  n  =  0  and 
n  =  which  giv 

iMa  X 

22  and  23,  respectively. 


o  =  n  -  which  gives  the  approximate  expressions  written  in  Equations 
m3  x 


dn  >  ^ 

dy‘ 


n=0 


-  1  i 

D  M 


(22) 


H—)2 

aldy' 


n=nmax  -  4An  "  ^  Ak}  "  nK,]  {23) 


The  derivative  in  Equation  22  can  be  evaluated  using  the  y  values  obtained 

previously.  The  derivative  in  Equation  23  can  be  specified  by  choosing  an 

nm  based  on  the  placement  of  the  outer  boundary.  Then  Equations  22 
max 

and  23  can  be  combined  to  give  the  following  equation  for  the  damping 
factor  D: 


D  =  ln{20 


{^) 

n=0  'dy; 


-2 


(24) 


With  the  damping  factor  D  known,  the  analytic  expression  (Equation  21) 
can  be  evaluated  on  each  interval  of  successive  y  values  obtained  from 
the  Blasius  solution  along  with  the  last  An  interval  near  the  outer  boun¬ 
dary.  This  procedure  gives  a  system  of  K  equations  for  the  K  unknowns 
Ak  which  are  then  solved.  Therefore,  values  for  the  amplification  factors 

and  damping  factor  are  obtained  at  each  e;  line  which  are  used  to  compute 
a  body  fitted  grid  system. 
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3.  BODY  SURFACE  POINT  DISTRIBUTION 

The  body  surface  can  be  defined  by  either  a  given  set  of  points  or 
by  an  analytical  expression.  The  body  surface  grid  points  can  then  be 
determined  by  either  interpolation  or  evaluation  of  the  analytical  expres 
sion. 


The  NACA  0012  airfoil  section  is  used  in  this  investigation.  This 
airfoil  has  been  used  in  several  investigations  (Reference  65)  and  is  an 
AGARD  (Reference  66)  designated  test  airfoil  section.  The  NACA  four  digit 
series  airfoils  have  an  analytic  description  of  the  body  surface  given  by 
(Reference  67) 

y  (x)  =  5t(0.2969x^  -  0.126x  -  0.3516x2  +  0.2843x3  -  0.1015x4}  (25) 

where  x  is  the  distance  along  the  chord  (nondimensional ized  by  chord),  yt 

is  the  upper  surface  ordinate,  and  t  is  the  thickness  given  as  a  fraction 
of  the  chord.  For  the  NACA  0012  airfoil,  t  =  0.12.  The  thickness  function 
given  by  Equation  25  does  not  close  the  body  at  the  trailing  edge  (x=l). 

A  circular  arc  is  used  to  close  the  body.  The  arc  is  constructed  to  have 
a  center  on  the  chord  line,  to  have  points  of  tangency  with  the  upper  and 
lower  airfoil  surfaces  near  x  =  0.99,  and  to  intersect  the  chord  line  at 
x  =  1 .0. 

The  distribution  of  grid  points  on  the  airfoil  surface  is  based  on 
both  resolving  streamwise  flow  field  gradients  and  defining  the  surface 
curvature.  These  criteria  require  the  clustering  of  points  in  the  leading 
edge  and  trailing  edge  regions  of  the  airfoil.  The  following  analytic 
transformation  is  used  to  obtain  the  desired  clustered  distribution  of  x 
abscissa  values 

tanh(z- ~^A-1- )  -  tanh(-  jg) 

tanh(1irL)  ■  tanh(- 

where  z  represents  the  input  set  of  equally  incremented  values  (0  z  1), 
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x  represents  the  clustered  set  of  abscissae  (0  <  x  <_  1),  and  A1 ,  A2,  and  A3 
are  arbitrary  constants.  The  constant  A1  determines  the  center  of  the 
hyperbolic  function  and  is  used  to  cluster  points  toward  one  end.  The 
constant  A2  varies  the  slope  of  the  function  at  the  center  and  thus 
determines  the  distribution  near  mid  chord.  The  exponent  A3  can  also 
affect  the  clustering  of  points  at  the  ends.  If  A3  >  1,  more  points 
occur  near  x  =  0;  while  if  A3  <  1,  points  are  clustered  more  at  x  =  1 . 
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SECTION  III 
GOVERNING  EQUATIONS 

The  time  dependent  incompressible  viscous  Navier-Stokes  equations 
for  two-dimensional  flows  are  presented.  The  primitive  variables  of  fluid 
velocity  and  pressure  are  used.  The  Reynolds  averaged  Navier-Stokes 
equations  are  obtained.  Boundary  and  initial  conditions  are  given. 

A  modified  turbulence  model  based  on  Prandtl's  mixing  length  theory  is 
described  which  expresses  the  Reynolds  stresses  in  terms  of  mean  flow 
quantities  for  separated  adverse  pressure  gradient  flows.  An  expression 
for  the  force  exerted  by  the  flow  field  on  the  body  is  derived.  The 
equations  and  boundary  conditions  are  transformed  by  the  curvilinear 
transformation  discussed  in  Section  II  and  written  in  the  computational 
plane. 

1 .  BASIC  EQUATIONS  OF  MOTION 

The  time  dependent  incompressible  Navier-Stokes  equations,  which 
describe  conservation  of  linear  momentum,  in  orthonormal  indicial  notation 
are  given  by 

8V 

+  3.(V.V.)}  =  -3.p  +  3.  t  (27) 

Jjl  1  J  J  ' 

where  p  is  the  fluid  density,  p  is  the  fluid  static  pressure,  V..  are  the 

components  of  fluid  velocity,  and  is  the  viscous  stress  tensor.  The 

subscript  i  =  1  corresponds  to  the  x  direction  and  i  =  2  to  the  y  direction. 
Conservation  of  mass  for  an  incompressible  flow  is  given  by 

3 .  V.  =  0  (28) 

J  J 

Turbulent  flow  is  introduced  by  assuming  that  the  flow  quantities  can 
be  described  in  terms  of  mean  and  fluctuating  quantities  in  the  form 
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where  the  bar  (-)  indicates  a  mean  quantity  and  the  prime  (')  indicates  a 
fluctuating  quantity.  In  this  form,  a  time  average  of  a  fluctuating 
quantity  is  assumed  to  be  zero.  Substitute  these  expressions  into 
Equations  27  and  28  and  time  average  each  equation  to  obtain  the 
Reynolds  averaged  Navier-Stokes  equations  for  two-dimensional 
incompressible  turbulent  flow  given  below. 


av. 

’{ar-  +3j(Vi)} 


-a.p  +  a  .{t..  -  p  v.  v.l 

v  j  ji  M  j  i 


(29) 


8  .  V.  =0  (30) 

J  J 

i — r 

The  term  -p  Vj  is  known  as  the  turbulent  Reynolds  stress.  The  incom¬ 
pressible  form  of  the  constitutive  equation  which  relates  the  stress  as 
a  function  of  the  rate  of  strain  is  now  introduced  as 


'ji  ■  “<sj  *  >i  V 


(31) 


where  y  is  the  dynamic  viscosity  coefficient.  The  turbulent  Reynolds 
stress  is  similarly  modelled  with  an  algebraic  eddy  viscosity  approach 
based  on  Prandtl's  mixing  length  theory.  The  corresponding  expression  is 

•  ■  '«,aj  7i  *  3i  V  1321 

where  eu  is  the  turbulent  eddy  viscosity.  In  this  way,  the  additional 

n 

— i - r 

unknown  -pV.  V.  is  expressed  in  terms  of  the  mean  flow  variables.  This 
J 

turbulence  closure  requires  an  expression  for  the  eddy  viscosity.  Models 
for  the  eddy  viscosity  are  examined  in  Section  III. 3. 


The  expressions  for  the  viscous  laminar  and  turbulent  stresses, 
(Equations  31  and  32)  are  substituted  into  Equation  29  to  obtain  the 


AFVIAL-TR-81  -3053 


simplified  turbulent  incompressible  Navier-Stokes  equations 


3V-  -  -  1  _  (v  +  eM) 

ttt1  +  a.(V  V  )  =  -  -  a  p  +  - 2- 

3t  J  J  1  pi  P 


'Vj  7f  4  sf  sj  V  (33) 


The  following  nondimensional  variables  are  introduced: 


x 


i 


P  = 


(34) 


where  c  is  the  airfoil  chord,  U  is  the  freestream  velocity,  and  p  is  the 

00  r  oo 

freestream  static  pressure.  Substituting  these  nondimensional  variables 
into  Equations  30  and  33,  with  the  time  average  and  nondimensional  symbols 
understood,  gives  the  following  set  of  equations 

D  =  3.  V  =  0  (35) 

J  J 


a  .(v.v.) 

j  j  i 


-3i 


p  + 


i 

Re, 


(a.a.Vi  +  3i3jV 


(36) 


where  D  is  the  divergence  of  velocity  and  Ret  is  the  turbulent  Reynolds 
number.  The  turbulent  Reynolds  number  is  defined  as 

Re  =  -ii -  (37) 

z  eM 

1  +  - 
M 


where  Re  =  pU^c/p,  the  Reynolds  number. 

Another  simplification  would  be  to  apply  Equation  35  into  Equation 
36.  The  divergence  D  is  contained  in  both  the  convective  and  viscous 
terms.  Although  analytically  correct,  the  velocity  divergence  is  in 
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general  not  identically  zero  in  a  numerical  result.  The  approach  used  by 
Hodge  (Reference  53)  and  developed  by  Harlow  and  Welch  (Reference  51)  is 
to  keep  the  divergence  D  in  the  convective  term  as  a  correction  term. 

This  approach  is  used  here  so  that  the  viscous  divergence  term  is  set 
equal  to  zero  in  Equation  36.  Since  the  convective  term  is  important 
throughout  the  '-'low  field  while  the  viscous  term  is  important  only  near 
the  airfoil  surface,  it  appears  reasonable  to  drop  the  viscous  divergence 
term  and  keep  the  convective  divergence  term  as  a  correction  term.  Thus, 
the  appropriate  set  of  equations  becomes 

sr-'VYi’  =  -3i  (38) 


D  = 


8  .V. 
J  J 


(39) 


The  calculation  of  pressure  in  the  primitive  variables  formulation 
must  be  carefully  handled.  The  pressure  does  not  appear  in  the  continuity 
equation  nor  does  a  time  derivative  of  pressure  occur.  Two  approaches 
have  been  extensively  used  to  overcome  this  difficulty.  Both  techniques 
compute  the  pressure  by  utilizing  the  velocity  divergence  D  as  a  correction 
factor.  The  calculated  pressure  is  made  to  satisfy  continuity.  The  first 
method,  introduced  by  Chorin  (Reference  68),  computes  pressure  using  an 
iteration  procedure  with  0  as  the  correcting  term.  The  expression  is 


Pls*’>  .  p<s>  .  ,0 


(40) 


where  s  denotes  the  iteration  number  and  $  is  an  acceleration  parameter 
which  can  be  either  a  constant  or  can  be  related  to  a  successive-over¬ 
relaxation  (SOR)  solution  of  a  Poisson  equation  for  pressure  as  shown  by 
Hodge  (Reference  53).  The  second  method  uses  a  Poisson  equation  for 
pressure  derived  from  the  divergence  of  the  momentum  equation  expressed 
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in  Equation  38.  The  form  of  the  equation  is  shown  in  Equation  41. 


D 


t 


=  7 


(41) 


where  a  represents  the  convective  and  viscous  terms.  If  the  expression 
is  algebraically  simplified  with  0  used  wherever  possible,  the  following 
simplified  incompressible  Poisson  pressure  equation  is  obtained 

Dt  =  -  {  (uD)x  +  (vD)^  +  u2  +  2vxuy  +  v2  +  uDx  +  vD^} 


~~  (D  +  D  }  -  (p  +  p  ) 
Re^  xx  yy  Kxx  Kyy 


(42) 


This  equation  contains  both  spatial  and  time  derivatives  of  the  velocity 
divergence  as  well  as  the  fluid  velocities  and  pressure.  Again,  continuity 
is  used  as  a  correcting  factor  for  calculating  pressure.  Further  simplifi¬ 
cation  can  occur  if  small  variations  in  the  velocity  divergence  are  assumed. 
In  this  case,  the  spatial  derivatives  of  D  are  set  to  zero  in  Equation  42 
and  the  result  becomes 

D.  +  (u2  +  2v  u  +  v2)  =  -  ' +  p  )  (43) 

t  '  x  x  y  y '  ..x  ryy' 

Hodge  (Reference  53)  used  both  Equations  42  and  43  and  found  no  apparent 
difference  in  the  flow  field  solutions  for  his  case  of  laminar  flow. 

Also,  Chorin’s  method  of  Equation  40  is  related  to  the  pressure  Equation 
43.  Hodge  (Reference  53)  showed  that  the  Chorin  technique  is  approximately 
related  to  a  solution  of  the  Poisson  pressure  equation  using  SOR  iteration 
on  a  coarse  grid  of  twice  the  spacing.  The  Poisson  pressure  equation 
technique  is  chosen  for  the  interior  flow  field  because  of  the  increased 
coupling  which  occurs  in  a  finite  difference  representation.  However, 
the  body  surface  pressure  is  calculated  using  the  Iteration  technique  of 
Equation  40.  Here  mass  conservation  is  Imposed  directly  and  the  difficulty 
of  evaluating  the  Laplacian  of  pressure  at  the  body  surface  in  Equation 
43  is  avoided. 
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The  following  set  of  equations  taken  from  Equations  38,  39,  and  43 
becomes  the  system  used  for  the  solution  of  incompressible  turbulent 
viscous  two-dimensional  flow  over  airfoils: 


u.  +  uu  +  vu  +  uD  =  -p 
t  x  y  Kx 

vt  +  uvx  +  vvy  +  vD  =  -p 
Dt  +  ux  +  2Vy  +  Vy  =  -(P 


+  Re^  ^uxx  +  uyy^ 

(44) 

+  Re^  ^vxx  +  vyy^ 

(45) 

xx  +  Pyy  ^ 

(46) 

(47) 

The  numerical  solution  of  the  governing  equations  is  obtained  by  per¬ 
forming  the  computations  in  the  computational  plane  obtained  with  the 
curvilinear  transformation  discussed  in  Section  II.  This  method  requires 
that  the  equations  are  transformed  to  the  computational  plane.  The 
governing  equations  (Equations  44  through  47)  are  transformed  using  the 
relations  given  in  Section  II  and  Appendix  A.  The  resulting  set  of  cor¬ 
responding  transformed  equations  becomes: 


*  0  <J,n  UC  ■  y(  un>  *  J  (x5  “n 


x  ur)  +  uD  = 

n  £ 


3  <yn  h  ■  h  p.’  *  rs: 


&s 


-  20U_  +  yU 

_Sn 


X 


nn) 


U  (-f) 

n  ,4 


+  u  (-4)} 


(48) 


vT  +  4  (y  v.  -  y  v  )  +  4  (x..  v  -  x  v.)  +  vD  = 
i  i  (av  -  20v_  +  yv  ) 

_  —  (  y  n  _  v  n)  +  — - —  { _ ^ 

J  vx£  pn  n  V  Re.  1  X^ 


+  %<X>  +  Vi)} 

J  J 


(49) 
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°T  *  7  t(»„  “t  ■  h  \)Z  *  2t,n  v5  '  yc  V  l*E  un  '  \  u5> 


,  ap  -  2gpr  +  yP 

*<*,  v  -  X„  »/>  .  -  (—a - -in - 21 

5  n  n  t,  jd 


+  Pn  [f  +  P?(7)} 


(50) 


D  "  (^n  UC  -  Un}  +  {\  \  ~  \  h] 


(51) 


2.  BOUNDARY  AND  INTIAL  CONDITIONS 

The  complete  definition  of  the  flow  problem  requires  the  specification 
of  sufficient  boundary  and  initial  conditions.  The  Navier-Stokes  equations 
are  coupled,  nonlinear,  first  order  in  time,  second  order  in  space,  elliptic 
partial  differential  equations.  The  elliptic  nature  of  the  equations  for 
velocity  requires  boundary  conditions  at  each  boundary.  The  first  order 
derivatives  of  pressure  require  a  boundary  condition  describing  the  free- 
stream  pressure.  The  time  derivative  terms  require  initial  conditions  on 
velocity  and  pressure  throughout  the  flow  field. 

Boundary  conditions  on  the  airfoil  surface  are  considered  rst. 

The  continuum  no  slip  boundary  condition  for  a  viscous  fluid  at  a 
stationary  solid  surface  is  imposed  on  the  tangential  velocity  at  the 
surface.  In  addition,  the  normal  component  of  velocity  must  be  zero 
for  a  stationary  surface  with  no  transpiration.  These  conditions  are 
met  by  the  specification  of  the  Dirichlet  boundary  conditions 

ufc  =  0  vb  =  °  (52) 

where  the  subscript  b  denotes  body  surface.  The  pressure  on  the  airfoil 
surface  is  computed,  and  no  body  boundary  condition  is  required  by  the  first 
derivative  pressure  terms. 
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The  far  field  boundary  conditions  for  velocity  and  pressure  specify 
the  freestream  flow  field.  Let  r  be  a  position  vector  in  the  flow  field 
with  the  origin  at  the  airfoil  mid-chord.  Then,  as  the  distance  |r|  from 
the  airfoil  becomes  large  the  fluid  velocity  and  pressure  approach  their 
nondimensional  freestream  values  as  follows: 


lim  (ui  +  vj)  =  cos  a  i  +  sin  a  j  (53) 


lim  p  =  0  (54) 


where  a  is  the  angle  of  attack  of  the  airfoil  chord  with  respect  to  the 
freestream. 

The  first  order  time  derivative  terms  require  initial  conditions  for 
the  velocity  and  pressure  in  the  field.  From  a  theoretical  viewpoint, 
since  the  asymptotic  steady  solution  is  desired,  any  continuous  initial 
specification  of  velocity  and  pressure  (subject  to  the  prescribed  boundary 
conditions)  is  permissible. 


The  boundary  conditions  and  initial  conditions  are  also  transformed 
to  the  computational  plane.  The  values  for  each  condition  remain  unaltered 
and  are  thus  not  repeated  here.  The  location  in  the  transformed  plane 
depends  on  the  specific  transformation  chosen  to  generate  the  computational 
grid.  This  matter  is  discussed  in  Section  IV  where  the  numerical  procedures 
are  described.  The  sensitivity  analysis  of  the  outer  boundary  location  is 
addressed  in  Appendix  F. 

3 .  TURBULENCE  MODEL 

Turbulence  modelling  is  of  considerable  interest  in  theoretical, 
experimental,  and  computational  research  work.  Many  aerodynamic  problems 
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of  practical  interest  contain  regions  of  turbulent  flow.  Several  models 
have  been  developed  which  vary  in  complexity  from  Prandtl's  mixing  length 
theory  to  techniques  employing  auxiliary  differential  equations  for  mixing 
lengths,  viscosity  coefficients,  and  turbulent  kinetic  energies.  The 
accuracy  of  each  model  varies  with  both  the  geometry  and  the  flow  field 
quantity  used  for  comparing  results.  In  computational  work,  the  ability 
to  translate  the  turbulence  model  into  numerical  procedures  and  the  computer 
resources  required  are  additional  factors  which  need  to  be  considered. 

In  this  investigation,  turbulence  is  modelled  with  an  algebraic  eddy 
viscosity  approach  based  on  Prandtl's  mixing  length  theory.  The  basic 

i  5 

concept  of  this  method  is  that  the  turbulent  Reynolds  stresses,  -pV.V., 

J  ' 

can  be  related  to  the  strain  rate  of  the  mean  flow  in  a  way  analogous  to 
that  for  laminar  viscous  stresses.  The  fundamental  expression  was  previously 
introduced  in  Equation  32  to  obtain  the  simplified  turbulent  incompressible 
Navier-Stokes  equation  (Equation  33).  The  algebraic  eddy  viscosity  has 
been  used  extensively  in  obtaining  numerical  solutions  for  a  variety  of 
aerodynamic  problems  (References  24,25,69). 

The  eddy  viscosity  model  for  turbulent  flow  with  zero  pressure  gradient 
near  a  solid  surface,  which  incorporates  the  modifications  of  Cebeci  and 
Smith  (Reference  63),  Van  Driest  (Reference  23),  and  Deiwert  (Reference 
57),  is  given  below  and  then  discussed. 

-2  -2  ^2 

e.  =  p(k1yDL)2  (u  +  vx)  (55) 


co  =  pk2ue5  Y 

0.  =  1  -  exp  {-  (^-M 

vA  p 


y  =  (1  +  5.5  (y/5)  } 


*  £ 
6  =  / 


/  0  -  *r)  dy 

Jo  u„ 
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The  tangential  surface  direction  is  given  by  x  and  the  normal  surface 
direction  by  y.  The  inner  and  outer  eddy  viscosities,  and  £q,  model 

the  Law  of  the  Wall  and  Law  of  the  Wake  (Reference  70)  regions.  The  inner 
mixing  length  is  assumed  proportional  to  the  distance  y  from  the  surface. 
k-|  is  the  Von  Karman  constant.  The  Van  Driest  (Reference  23)  damping  factor 

models  the  laminar  sublayer  region  where  is  the  surface  shear  stress  and 

+  * 

A  is  the  decay  constant.  The  displacement  thickness  6  defined  by  Equation 

59  is  the  outer  layer  length  scale.  The  intermittency  factor  y»  defined 

in  Equation  58,  models  the  decrease  of  the  turbulence  intensity  from  a 

fully  turbulent  region  to  the  inviscid  region  away  from  the  body  surface. 

The  intermittency  is  a  function  of  the  normal  distance  from  the  wall 

nondimensional ized  by  the  local  boundary- layer  distance  <5.  The  local 

edge  boundary-layer  velocity  ug  is  the  outer  layer  characteristic  velocity. 

The  standard  Cebeci-Smith  (Reference  63)  values  for  the  constants  are  k^ 

=  .41,  =  .0168,  and  A+  =  26.  Continuity  of  the  eddy  viscosity  distribu¬ 

tion  is  achieved  by  switching  from  the  inner  to  the  outer  model  whenever 

e.  first  exceeds  e  . 
l  o 

Recently,  Jobe  and  Hankey  (Reference  71)  compared  this  eddy  viscosity 
model  given  by  Equations  55  through  59  with  experimental  data  for  several 
constant  adverse  pressure  gradient  attached  flows.  They  varied  the  k^ , 

k£,  and  A+  parameters  in  a  numerical  turbulent  boundary-layer  program  until 

the  solutions  matched  the  experimental  boundary-layer  velocity  profiles 
and  skin  friction  coefficients.  These  parameters  had  substantially  different 
values  than  the  previously  cited  values  obtained  for  zero  pressure  gradient 
flows.  They  found  that  when  the  flow  encounters  an  adverse  pressure  gradient 
the  inner  eddy  viscosity  parameters  rapidly  adjust  to  the  values  of  k-j  = 

.65  and  A+  =  52.  The  outer  eddy  viscosity  constant  appears  to  linearly 
increase  with  distance  measured  from  the  start  of  the  adverse  pressure 
gradient.  Thus,  Jobe  and  Hankey  (Reference  71)  found  that  k^  can  vary 

from  a  value  of  .0168  to  values  of  .03  or  greater  given  that  the  adverse 
pressure  gradient  exists  for  a  sufficient  length  downstream. 


31 


AFWAL-TR-81 -3053 


The  use  of  an  eddy  viscosity  model  similar  in  form  to  the  one  given 
in  Equations  55  through  59  for  an  airfoil  at  angle  of  attack  requires 
examination  of  the  typical  flow  field  characteristics.  The  upper  surface 
of  the  airfoil  has  a  pressure  distribution  which  includes  a  region  of  low 
pressure,  the  suction  peak,  near  the  airfoil  leading  edge.  The  suction  peak 
is  followed  by  an  adverse  pressure  gradient  region.  The  pressure  gradient 
tends  to  decrease  as  the  flow  recovers  toward  the  trailing  edge.  The 
turbulence  model  that  is  used  for  an  upper  surface  turbulent  boundary-layer 
should  include  the  effects  of  this  adverse  pressure  gradient.  The  inner 
eddy  viscosity  parameters  take  on  the  Jobe  and  Hankey  (Reference  71) 
values  of  k^  =  .65  and  A+  =  52  in  this  region.  The  outer  eddy  viscosity 
parameter  k^  may  also  change  in  the  adverse  pressure  gradient  region. 

However,  the  rate  of  change  is  affected  by  the  pressure  gradient  and 
downstream  distance.  The  pressure  gradient  does  decrease  rapidly  as  the 
flow  proceeds  downstream.  Any  tendency  for  k2  to  increase  is  probably 

offset  by  the  decreasing  adverse  pressure  gradient.  Thus,  the  value  of 

is  kept  at  the  equilibrium  value  of  .0168. 

The  eddy  viscosity  model  with  the  discussed  adverse  pressure  gradient 
modifications  is  now  written  in  terms  of  the  nondimensional  variables 
introduced  in  Equations  34  and  37.  All  length  quantities  are  nondimen- 
sionalized  with  respect  to  the  airfoil  chord.  The  mean  and  nondimensional 
quantities  are  henceforth  implied.  The  eddy  viscosity  model  becomes 


T  *  <  W2  <uy  +  Re 

T  =  kzV*  Re 

\  •  1  -  exp  {-  (|uyb|Re)  } 
A 

-1 

y  -  (1  +  5.5  (y/6)6} 


(60) 

(61) 

(62) 

(63) 
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R 


(64) 


where  ,  k^. 


and  A 


are  considered  adverse  pressure  gradient  parameters. 


The  eddy  viscosity  model  must  also  represent  the  physical  characteristics 
in  large  separated  regions  which  occur  during  airfoil  stall.  Trailing  edge 
stall  occurs  when  the  turbulent  boundary-layer  separation  moves  toward  the 
leading  edge  of  the  airfoil  as  the  angle  of  attack  increases.  The  pressure 
distribution  in  a  trailing  edge  separated  region  has  small  gradients 
(Reference  72).  This  observation  motivates  another  change  to  the  eddy 
viscosity  model.  The  inner  eddy  viscosity  parameter  is  relaxed  back 
toward  the  equilibrium  value  of  .41.  The  separation  point,  which  precedes 
the  return  to  zero  pressure  gradient,  is  selected  as  the  point  to  start 
the  exponential  decay.  The  following  relaxation  form  with  distances 
nondimensional ized  by  the  chord  is  used 


f  =  A{1  +  B  exp( - - — -)  s  -  sd 

f  =  1  s  <  sd 


(65) 


where  f  is  the  relaxation  factor  which  multiplies  the  initial  value  of 
2  2 

k-j ,  k^,  s  is  the  distance  downstream  from  the  trailing  edge  separation 

point,  Sj  is  the  delay  distance,  sf  is  the  relaxation  distance,  and  A  and 

B  are  constants.  The  constants  A  and  B  are  determined  by  applying  the 

2  2 

conditions  f  =  1  at  s  =  sd  and  f  =  k^/k^  as  s  approaches  infinity  where 
k^  is  the  asymptotic  value  of  k^  downstream.  Thus,  A  and  B  become 

A  =  k^f/k*.  B  =  A’1  -1  (66) 


a 
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Measurements  of  Bachalo  (Reference  73)  and  Baker  (Reference  74)  in 
the  trailing  edge  separation  region  indicate  the  presence  of  a  lower  but 
finite  turbulent  intensity  when  compared  with  the  separated  turbulent 
shear  layer.  The  velocity  derivative  terms  in  the  inner  eddy  viscosity 
and  damping  factor  expressions  (Equations  60  and  62)  can  approach  zero 
in  a  region  near  a  velocity  inflection  point.  In  this  case,  a  zero  eddy 
viscosity  is  permitted  in  a  turbulent  separated  region  which  is  contrary 
to  the  apparent  flow  structure.  This  same  zero  eddy  viscosity  behavior 
can  also  occur  in  the  case  of  shedding  separation  bubbles  from  the  airfoil 
surface.  In  order  to  prevent  the  appearance  of  laminar  regions  within 
the  turbulent,  adverse  pressure  gradient,  separated  flow  field  during 
airfoil  stall,  a  limiting  technique  is  employed.  The  inner  eddy  viscosity 
is  prevented  from  decreasing  in  value  as  the  normal  distance  from  the 
surface  increases.  This  restriction  simulates  the  uniform  turbulent 
intensities  measured  (References  73,74)  in  the  trailing  edge  separated 
region.  The  inner  eddy  viscosity  is  further  limited  by  imposing  a  condition 
of  no  decrease  in  the  downstream  direction.  This  limit  provides  for  a 
finite  turbulent  intensity  in  separated  regions  which  may  develop  in  the 
turbulent  boundary  layer.  The  limit  is  initiated  by  obtaining  the  distri¬ 
bution  of  eddy  viscosity  in  the  attached  boundary-layer  near  the  leading 
edge  of  the  airfoil.  This  distribution  is  compared  with  the  distri¬ 
butions  calculated  by  the  model  at  downstream  locations.  The  outer  eddy 
viscosity,  which  has  no  velocity  derivative,  is  not  limited  in  any  way. 

The  "limiting"  technique  is  passive  in  the  sense  that  the  locally  computed 
value  for  the  inner  eddy  viscosity  is  used  whenever  possible.  Numerical 
experiments  without  the  "limiting"  technique  displayed  unphysical  results 
because  the  conventionally  modelled  eddy  viscosity  is  divergent.  The 
inception  of  separated  flow  decreased  the  eddy  viscosity  thereby  causing 
futher  separation.  The  "limiting"  concept  is  necessary  to  prevent  this 
unphysical  result. 
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The  transition  from  laminar  flow  to  fully  turbulent  flow  is  modelled 
using  the  transition  model  of  Dhawan  and  Narasimha  (Reference  75).  The 
expression  for  the  transition  factor  r  is  given  below 

x  -  x  2 

r  =  1  -  exp{-  . 41 2 ( - }  (67) 

A 

where  x  is  the  tangential  direction  coordinate  on  the  body  surface,  x^.  is 
the  location  where  transition  starts,  and  A  =  xr=3^4  -  xr=i/4  defines  the 
distance  required  to  proceed  from  r  =  1/4  to  r  =  3/4. 

Some  of  the  expressions  in  the  turbulent  boundary-layer  eddy  viscosity 
model  are  conveniently  expressed  in  terms  of  quantities  in  the  transformed 
plane  described  in  Section  II.  The  velocity  derivatives  in  the  inner  eddy 
viscosity  model  are  formulated  by  determining  the  tangential  and  normal 
components  of  velocity,  u  and  v  respectively,  to  the  body  surface  n  =  n-j 

at  a  position  £ 


u(e,  n)  =  c2u(c,  n)  -  c.jv(s,  n) 

(68) 

v(C,  n)  =  c]u(c,  n)  +  c2v(s,  n) 

(69) 

where  c^  and  c2  are  components  of  the  unit  outward  normal  to  the  body 
surface  at  position  {?,  n-| )  given  by 

Cj  =  -y  c2  =  x  ^//y 

Then,  the  directional  derivative  definition  is  used  to  obtain  the  surface 
normal  derivative  of  u  and  the  surface  tangential  derivative  of  v 
expressed  as 

Is-  =  {cl(Ve  '  +  C2(Vn  '  V?)}/J  (70) 


35 


AFWAL-TR-81 -3053 


31 


f?  *  (c2(Vs  -  yt%>  '  cl(Yn  ‘  V<L),/J  (7'> 

which  correspond  to  the  and  derivatives  of  the  Cartesian  eddy 
viscosity  model  expressions  {Equations  60  through  64). 

The  turbulence  structure  in  the  far  wake  is  modelled  with  an  eddy 
viscosity  expression  similar  to  the  outer  eddy  viscosity  model  in  the 
turbulent  boundary-layer.  The  model  is  based  on  the  assumption  of  a 
sel f- preserving,  equilibrium  turbulent  wake  with  constant  pressure  in 
incompressible  flow.  The  approach  has  been  experimentally  investigated 
by  Narasimha  and  Prabhu  (Reference  76)  who  obtained  the  following  expression 

e  =  pk,w  6  (72) 

w  3  o  w 

where  wq(x)  =  max  {U^  -  u(x,y)}  is  the  maximum  velocity  defect  in  the  wake, 

6  is  a  half  wake  width  defined  where  the  velocity  defect  becomes  one  half 
w 

the  maximum  value,  and  k^  is  a  constant  equal  to  .065.  If  the  experimentally 

integrated  value  of  the  wake  defect  function  is  used,  then  Equation  72  can 
be  written  in  terms  of  the  parameters  of  edge  velocity  and  displacement 
thickness  similar  to  boundary-layers  as  follows: 

ew  =  k,u  $  Re  (73) 

—  3  e  w 

u 

* 

where  u  is  the  nondimensional  velocity  at  the  edge  of  the  wake,  6  is  the 

nondimensional  displacement  thickness  for  the  half  wake  as  defined  by 
Equation  64,  and  k^  becomes  .0634.  This  expression  has  been  utilized  by 

Green,  et  al  (Reference  77)  in  an  integral  method  for  predicting  two-dimen¬ 
sional  Incompressible  and  compressible  turbulent  wakes  of  lifting  airfoils. 
For  an  incompressible,  constant  pressure  wake  the  momentum  integral  equation 
reduces  to  a  statement  that  the  momentum  thickness  is  a  constant.  These 
conditions  exist  in  the  far  wake  of  an  airfoil.  In  addition,  in  the  far 
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wake  the  displacement  and  momentum  thicknesses  become  equal  as  shown  by 
the  experimental  results  of  Narasimha  and  Prabhu  (Reference  76).  Thus 
the  far  wake  eddy  viscosity  is  approximately  constant.  The  transition 
in  the  wake  from  the  turbulent  core  region  to  the  inviscid  freestream 
region  on  each  side  of  the  wake  is  modelled  with  the  intermittency  factor 
of  Equation  63  where  <5  is  now  the  wake  half  width  and  y  is  the  normal 
distance  from  the  wake  center. 

The  turbulent  near  wake  region  mixing  rates  where  the  wall  boundary- 
layers  merge  and  eventually  become  the  far  wake  are  not  well  understood 
from  either  an  analytical  or  experimental  viewpoint.  In  the  flow  regime 
under  consideration,  a  laminar  boundary-layer  from  the  lower  airfoil  surface 
mixes  with  a  turbulent  boundary-layer  from  the  upper  surface.  An  eddy 
viscosity  model  similar  to  that  used  in  the  outer  turbulent  boundary-layer 
and  the  far  wake  is  used.  Inouye,  et  al  (Reference  78)  applied  the 
Cebeci-Smith  (Reference  63)  outer  layer  model  in  the  near  wake  of  a  flat 
plate  and  obtained  good  agreement  for  the  mean  velocity  profiles.  They 
allowed  the  constant  k^  to  reach  twice  the  equilibrium  value  of  .0168 

over  a  length  of  one  chord  from  the  trailing  edge.  The  computed  eddy 
viscosity  was  still  about  50%  low  when  compared  with  experiment  at  the 
one  chord  distance.  Later,  Burggraf  (Reference  79)  compared  several 
turbulence  models  for  various  near  wake  flows  and  concluded  that  the 
Cebeci-Smith  model  adequately  predicted  the  mean  velocity  profiles. 


The  merging  of  the  upper  and  lower  surface  boundary-layers  is  modelled 
with  a  simplified  interaction  hypothesis  approach  introduced  by  Bradshaw 
(Reference  80).  The  two  layers  are  initially  allowed  to  develop  downstream 
a  distance  WL^  without  interaction.  At  this  point  the  eddy  viscosity  is 

increased  exponentially  in  the  adjacent  laminar  boundary-layer  from  a 
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zero  value  to  the  far  wake  calculated  value  e  at  a  distance  WL0  from 

ow  2 

the  trailing  edge.  This  eddy  viscosity  variation  is  formulated  as  follows 

s  -  WL 

e  =  1/2  eow  {1  +  tanh  (8  WL - wf  )}  s  >  WL1 

ow  wl2  wl1  1  (74) 

e  =  0 


s  <_  WL1 


where  WLay  is  the  distance  from  the  trailing  edge  to  a  location  half  way 

between  the  distances  WL-|  and  WL2  and  s  is  the  local  distance  from  the 

trailing  edge.  The  coefficient  8  is  selected  so  that  e  =  0  at  s  =  WL^ 

and  e  =  eQw  at  s  =  WL2  because  at  s  =  WL2  the  hyperbolic  tangent  function 

becomes  equal  to  0.9993  while  at  s  =  WL^  it  becomes  -  0.9993.  The  eddy 

viscosity  profile  in  the  turbulent  boundary-layer  near  the  trailing  edge 
is  extended  into  the  near  wake. 


4.  FORCE  COEFFICIENTS 

The  force  which  is  exerted  on  the  body  by  the  moving  fluid  around  the 
body  can  be  determined  from  the  Cauchy  integral  equation  for  conservation 
of  linear  momentum  applied  at  the  body  surface  or  on  a  closed  path  surround¬ 
ing  the  body.  The  basic  expressions  for  the  force  coefficients  formulated 
in  the  (x,y)  physical  plane  coordinate  system  are  given  by  Equations  D-12a 
and  D-12b  derived  in  Appendix  0  for  an  incompressible  turbulent  viscous  flow. 

Consider  first  the  body  surface  where  the  components  of  velocity  are 
identically  zero  from  the  no  slip  boundary  conditions.  In  this  case,  the 
equations  become 

Cfx  ■  /  2v  *  k  (4"1  £  *  2"2  <!y  *  *  |75> 

^b 

Cfy  ■  /  +  Ae  <2"1  *  fj>  ♦  4"2  4s  <76> 

where  the  integration  is  over  the  body  contour  denoted  by  S^.  The  first 
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term  in  each  expression  is  the  contribution  caused  by  pressure  forces  while 
the  second  term  is  the  viscous  stress  contribution. 


The  force  coefficients  (Equations  75  and  76)  determined  at  the  body 
surface  can  be  transformed  to  the  computational  plane  by  using  the  relations 
for  the  outward  normal  vector  to  an  n  contour,  line  integrals  and  derivatives 
given  in  Appendix  A.  The  transformed  equations  become 


/max  , 

{V  +  m  K  +  y<-(Vn  •  Vn))}dC  (77) 


^min 


C,  =  2 

fy 


/max  , 

‘-V  +  S5 

Sin 


ReJ  K  -  x5(Vn  '  Vn))}dC  (78) 


where  the  no  slip  boundary  condition  which  implies  =  v?  =  0  on  the  body 

surface  is  used.  The  viscous  stress  terms  are  further  simplified  by  using 
continuity  at  the  body  surface  and  the  vorticity  (_7  x  v)  Equation  A. 14  on 
the  body  surface  to  obtain 

/^max  2x  u> 

<  V  •  K*-  1  «  <79> 

S,fn 

imax  Zy  u 

Cf,  '/  '-V-RS5-  >  dE  (80) 

^min 


where  the  vorticity  on  the  body  surface  u  is  given  by 

i 

“  B  +  Vn)/J  (si) 
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The  corresponding  moment  (positive  counterclockwise)  about  a  point  P 
located  within  the  body  cross-section  with  coordinates  (Xp,  y  )  is  given 
as 

^max  2 y  u  2x  u 

CMP  =  J  <(x  -  xp)(-2x?p  -  -  (y  -  yp)(2ycp  -  — )}d£  (82) 

^min 


Next  consider  again  the  more  general  case  of  a  closed  path  around 

the  body  whose  force  coefficients  are  given  by  Equations  D-12a  and  D-12b. 

Assume  that  the  closed  path  is  an  n  contour.  Then,  using  the  transformed 

expressions  in  Appendix  A  for  the  outward  normal  vector  to  an  n  contour, 

integrals,  and  derivatives,  the  force  coefficients  for  the  body  in  terms 

of  quantities  on  a  constant  n  line  in  the  flow  field  become 

P 


max 


:fx  \ln  '2V  *  ^7  [U"(*E2  +  2y2)  '  Vxt\  *  +  vE(Vn> 


-V„(V5)]  *  2"!uyt  'VXE)1  dt  '  It 


(83) 


'fy 


2iriax 

J  {‘2xeP  +  Rrj  [vn(2xr  + 


^min 


■  Re  j  ‘  v^2x5xn  +  y£yn)  +  ue^xqy^ 


Ti„  S 


-  ^(x^)]  +  2v(uy?  -  vx5)>  dc 


/7max 

/ 


2v|j|d?dn  (84) 


nb  ^min 


where  nk<  n  <  n _ .  The  terms  in  the  line  integral  represent  in  order 

b  —  p  —  max  J 

pressure  forces,  viscous  forces,  and  convective  outflows  of  linear  momentum. 

The  area  integral  term  is  the  time  rate  of  increase  of  linear  momentum 

contained  in  the  control  volume  bounded  by  the  body  contour  and  the 

selected  contour. 

The  corresponding  lift  and  drag  coefficients  are  then  calculated  using 
Equations  D-13  and  0-14. 
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SECTION  IV 
NUMERICAL  METHOD 

The  numerical  methods  used  to  obtain  solutions  for  incompressible 
turbulent  flow  over  airfoils  are  presented.  The  numerical  finite  difference 
procedure  which  yields  the  grid  transformation  is  first  discussed.  The 
implicit  finite  difference  method  for  obtaining  approximate  solutions  to 
the  unsteady,  incompressible,  turbulent  Navier-Stokes  equations  is  next 
described.  Numerical  boundary  conditions  and  initial  conditions  are  con¬ 
sidered.  The  numerical  procedure  for  incorporating  the  turbulent  eddy 
viscosity  model  is  discussed.  The  computation  of  the  force  coefficients 
on  the  body  is  then  presented. 

1 .  GRID  TRANSFORMATION 

The  numerical  solution  of  the  governing  equations  for  the  transforma¬ 
tion  is  carried  out  in  the  transformed  or  computational  plane  on  a  square 
mesh  (5,  n)  system.  The  body  contour  Cb  and  the  far  field  boundary  C^. 

become  constant  n  contours  (Figure  1).  The  far  field  boundary  is  usually 
used  to  approximate  infinity  in  incompressible  flow  problems.  Constant  n 
contours  in  the  physical  plane  are  required  to  transform  to  equi-spaced  ri 
coordinate  lines  with  unit  step  size  in  the  transformed  plane.  The  spacing 
of  n  contours  in  the  physical  plare  is  controlled  by  the  set  of  elliptic 
equations.  The  £  contour  spacing  on  the  body  surface  and  outer  boundary 
in  the  physical  plane  is  determined  by  the  chosen  point  distribution  on 
each  surface.  The  constant  £  contours  are  also  required  to  transform  to 
equi-spaced  £  coordinate  lines  with  unit  step  size  in  the  transformed 
plane.  The  £  lines  are  denoted  by  the  subscript  i  with  range  1  i  i  f.  IMAX 
where  IMAX  is  the  number  of  £  lines.  Similarly,  the  n  lines  are  identified 
by  a  subscript  j  with  range  of  1  ±  j  ±  JMAX  where  JMAX  is  the  number  of  n 
lines.  The  values  j  =  1  and  j  =  JMAX  are  the  body  and  outer  boundary 
surfaces,  respectively.  The  constant  £  lines  i  =  1  and  i  =  IMAX  denote 
the  branch  cut  in  the  physical  plane  and  are  therefore  identical. 
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The  governing  elliptic  transformation  equations  (Equations  9a  and  9b) 
are  written  with  finite  differences  in  quasi-linearized  form  for  the 
location  (i,  j).  Second  order  accurate  central  differences  given  in 
Appendix  B  are  used  for  evaluating  each  term.  The  highest  order  derivative 
contains  the  unknown  at  location  (i,  j)  while  the  transformation  coefficients, 
which  contain  lower  order  derivatives,  are  evaluated  with  previous  values 
of  the  variables  x  and  y.  The  difference  expressions  are  given  below  where 
the  indices  i  and  j  are  understood  when  omitted. 


xi j  =  {a(xi+l  +  xi-l)  ‘  2  (xi+l  j+1  "  xi+l  j-1  +  xi-l  j-1  "  xi-l  j+1 J 


(85) 


+  +  )  *  2  C(xi+i  ”  x-j_i  )P  *  ~  **"  y) 


yi j  ^yi+l  +  yi-l^  “  2  ^yi+l  j+1  ”  yi+l  j-1  +  yi-l  j-1  '  yi -1  j+1  ^ 
+  Y(yj+1  +  yj-l5  +  T  C(yi+1  ‘  yi-l)P  +  (yJ+l  *  yj-l)Q]}/2(a  +  y) 


(86) 


The  i  index  near  the  branch  cut  is  adjusted  as  follows: 

when  i  =  1  then  i  -  1  =  IMAX  -  1 

when  i  =  IMAX  -  1  then  1+1=1 

when  i  =  IMAX  then  i  =  1 

The  set  of  finite  difference  equations  has  the  range  of  indices  given  by 

1  <  1  <  IMAX  -  1  and  2  <  j  <  JMAX  -1  which  results  in  (IMAX-1)  (JMAX-2) 
equations  for  the  unknown  values  of  x  and  y  at  the  interior  grid  points. 
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The  system  of  difference  equations  is  solved  using  successive  over 
relaxation  (SOR)  iteration.  For  a  general  function  f,  an  approximate 
solution  for  f  using  SOR  iteration  has  the  form 

f(s+1)  =  uf*  +  (1  -  w)f^  (87) 

where  s  represents  the  iterate  number,  u>  is  an  acceleration  parameter, 
and  *  represents  the  current  value  of  f  calculated  from  the  difference 
equation.  The  acceleration  parameter  is  always  one  for  the  first  iteration. 
The  system  of  linearized  difference  equations  given  by  Equations  85  and  86 
obtained  from  the  elliptic  generating  differential  equations  is  consistently 
ordered  (Reference  81).  Thus,  local  optimum  acceleration  parameters 
can  be  calculated  from  matrix  theory  (Reference  81)  and  are  given  by  Hodge 
(Reference  53).  Convergence  of  the  iteration  procedure  is  determined  by 
comparing  the  absolute  value  of  the  difference  in  successive  iterates 
with  a  prescribed  error  criteria  E. 

The  SOR  iteration  technique  requires  an  initial  guess  for  the  solution. 
Geometric  contours  with  a  shape  similar  to  the  outer  computational  boundary 
are  used  for  constant  n  lines  and  straight  lines  which  emanate  from  the 
body  center  are  used  for  ?  lines.  A  similar  approach  has  been  success¬ 
fully  used  by  Thames  (Reference  44)  and  Hodge  (Reference  53). 

The  forcing  function  Q(?,  n)  defined  by  Equation  15  is  evaluated 
throughout  the  U,  n)  plane  following  the  procedure  described  in  Section 
II. 2.  The  similarity  parameter  values  c  (which  give  equal  increments  of 
tangential  velocity  for  a  specified  number  of  boundary-layer  points) 
are  calculated  from  Equation  18  using  Newton-Raphson  iteration.  The 
equivalent  y  values  at  a  given  x(s  line)  are  found  from  the  similarity 
relation  (Equation  19).  The  n  derivatives  given  by  Equations  22  and  23 
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are  then  approximated  using  the  above  values  for  y  and  the  specified  n 
differences  as  follows: 


where  An  =  1 ,  y^  are  the  boundary-layer  y  values,  and  Ay^  is  the  specified 

far  field  y  difference  between  the  last  two  points  on  the  given  5  line. 

Then  the  damping  factor  D  is  found  using  Equation  24.  Next,  the  analytical 
expression  for  dn/dy  in  Equation  21  is  evaluated  on  the  far  field  interval 
A y^  and  on  each  boundary-layer  Ay  interval  determined  by  the  y^  values. 

This  system  of  equations  for  A ^(c)  is  solved  using  Gauss  elimination. 

Thus,  for  a  given  i  line,  values  for  Ak  and  D  are  calculated  for  use  in 
evaluating  the  Q(£,  n)  forcing  function.  This  function  is  then  used 
during  the  SOR  iteration  for  the  transformation  values  of  x  and  y. 

Additional  ?  contours  can  be  added  in  the  wake  region  by  using  the 
interpolating  technique  given  in  Appendix  C. 

2.  NAVIER-STOKES  FINITE  DIFFERENCE  EQUATIONS 

Approximate  solutions  of  the  unsteady  (Reynolds  averaged,  incompres¬ 
sible,  turbulent,  and  Navier-Stokes}  equations  are  obtained  by  using  an 
implicit  finite  difference  procedure.  Explicit  numerical  methods  can  be 
used  to  obtain  approximate  solutions  by  advancing  the  time  by  small 
increments.  The  allowable  time  increments  are  constrained  by  numerical 
stability  criteria.  Implicit  finite  difference  methods  can  also  provide 
approximate  solutions  by  marching  in  time  with  small  time  steps.  Implicit 
methods  usually  do  not  have  numerical  stability  constraints  imposed  on  the 
time  step.  However,  at  each  time  step  a  large  system  of  simultaneous 
difference  equations  must  be  solved.  A  matrix  iterative  technique  such 
as  SOR  iteration  is  an  attractive  solution  method  because  of  the  large 
matrix  size.  Convergence  criteria  are  thus  introduced  in  place  of 
stabil ity  criteria . 
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The  implicit  method  consists  of  both  a  linearization  of  the  equations 
and  the  application  of  finite  differences  to  those  equations.  The  x  com¬ 
ponent,  transformed,  momentum  equation  given  by  Equation  48  is  considered 
first.  Rearranging  the  convective  terms,  the  equation  becomes 


"T  *  ut(u  f  -  v  r>  *  un(v  f  -  "  J1’  +  uD  '  -hf  *  "n  / 

1  ,  '  2eUen  * 


+  Ret  <■ 


ja  *  U,  (4)  +  u  (-4)) 

?  Jz  n  Jz 


(89) 


Equation  89  is  now  written  in  quasi-! inear  form  where  lower  order 
derivatives  of  the  variables  occur  in  "coefficients".  The  current  value 
for  the  unknown  variable  at  location  (i,  j)  occurs  in  the  highest  order 
derivative  of  each  term  while  the  "coefficients"  are  evaluated  with  the  last 
available  values  for  the  variables  u,  v,  and  p.  In  this  form,  the  coefficient 
of  u^  in  the  convective  terms, j(uy^-vx^)/J,  is  proportional  to  the  component 

of  velocity  in  the  £  direction.  Similarly,  the  coefficient  of  u^  in  the 
convective  terms, ^(vx^-uy^)/J,  is  proportional  to  the  component  of  velocity 

in  the  n  direction.  These  quantities  are  designated  UC  and  VC,  respectively. 

Finite  difference  expressions  given  in  Appendix  B  are  used  to  approxi¬ 
mate  the  terms  in  Equation  89.  First  order  accurate  backward  differences 
are  used  for  the  time  derivative.  The  first  derivatives  of  pressure  are 
approximated  with  second  order  accurate  central  differences.  Likewise, 
second  order  central  difference  expressions  are  used  for  the  second 
derivatives  of  velocity  in  the  viscous  terms.  Second  order  accurate 
upwind  differencing  is  used  for  the  velocity  first  derivatives  in  the 
convective  terms.  The  quantities  UC  and  VC  determine  the  sign  of  the 
local  velocity  in  the  £  and  n  directions.  If  UC  >_  0,  then  backward 
upwind  differencing  is  used  for  the  £  derivatives  while  UC  <  0  dictates 
forward  upwind  differencing.  Also,  backward  upwind  differencing  for  n 
derivatives  occurs  if  VC  >_  0  while  forward  upwind  differencing  is  used 
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if  VC  <  0.  Second  order  accurate  upwind  differences  are  also  used  for  the 
velocity  first  derivatives  in  the  viscous  terms.  The  direction  is  determined 
by  the  sign  of  the  coefficients.  Points  on  the  J  =  2  and  J  =  JMAX  -  1 
contours  use  first  order  upwind  differences  for  n  first  derivatives  if 
the  velocity  direction  indicates  flow  from  the  body  surface  (J=l)  or  the 
outer  boundary  (J  =  JMAX).  Near  the  body  surface  the  grid  spacing  is 
extremely  small,  and  velocity  gradients  are  small  near  the  outer  boundary. 
Thus,  in  both  cases  the  associated  artificial  viscosity  has  a  negligible 
effect.  The  transformation  derivatives  in  Equation  89  are  evaluated  with 
second  order  central  differences.  The  velocity  derivatives  in  the  velocity 
divergence  term  are  also  evaluated  using  central  differences. 


Difference  equations  are  written  with  the  following  convenient 
conventions.  Terms  which  contain  unknowns  at  location  (i,j)  are  fully 
subscripted.  The  finite  difference  expression  for  location  (i,  j)  as 
described  previously  is  understood  when  only  the  term  is  given.  The 
subscripts  (i,  j)  and  superscript  n  are  assumed  when  omitted.  The 
resulting  difference  equation  for  Equation  89  becomes 


ui j  '  Ui j  +  (3u .  -  4uTr,  +  uTr7)|UC|  +  (3u . 


4uJCl  +  uJC2)1VCl+  uD 


=  '  23  (pi+l  "  pi-l*  +  2J  (pj+l  ‘  Pj-1 5 


+  *k{7 (Ui+1 '  2uij  +  ui-i}  +>{Vi  -  2uij  +  UJ-1} 


(ui+l  j+1  *  Vl  j+1  +  ui-l  j-1  '  Ui+1  j-l} 


(3ui  -  4uiyl  +  uIV2)jUV|  -  (3u.  -  4uJV]  +  uJV2)|VV|} 
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where 


for 

UC  >  0 

for 

UC  < 

:  0 

IC1 

=  i  -  1 

IC1 

=  i 

+  1 

IC2 

=  i-2 

IC2 

=  i 

+  2 

for 

VC  >  0 

for 

VC  « 

:  0 

JC1 

=  j  -  1 

JC1 

=  j 

+  1 

JC2 

=  j  -  2 

JC2 

=  J 

+  2 

and  UV  =  -t/2J2  and  VV  =  -a/2J2 
where 


for 

UV  >  0 

for 

UV  < 

:  0 

IV1 

=  i  -  1 

I VI 

=  i 

+  1 

IV2 

=  i  -  2 

IV2 

=  i 

+  2 

for 

VV  >  0 

for 

VV  « 

:  0 

JV1 

=  j  -  1 

JV1 

=  J 

+  1 

JV2 

=  j  -  2 

JV2 

=  J 

+  2 

The  difference  equation  for  conservation  of  mass  (Equation  51)  for  location 
(i,  j)  is  given  by 


D  =  2J  (ui+l  '  ui-1)  "  2J  (vi+l  '  vi-l)  "  2J  (uj+1  ‘  uj-l) 


(91) 


The  following  finite  difference  definitions  are  introduced  for  the 
transformation  derivatives  and  coefficients  at  location  (i,  j): 

XETA  =  x  /2J  XXI  =  x„/2j 

n  K 

YETA  =  yn/2J  YXI  =  y^/20 

ALFA  =  a/J2  GAMA  =  y/J2  BETA  =  -S/2J2 
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With  these  definitions.  Equation  90  can  be  written  to  provide  an  expression 

for  u".  given  by 
1  j  a 

Uij  =  {uij1  +  At[-YETA(pi+1  -  p1_1)  +  YXI(pj+1  -  p.^) 

+  |UC|(4uIcl  -  uIC2)  +  | VC | (4u JC1  -  uJC2)  -  uD 
+  Re^  [ALFA(ui+1  +  u. _1 )  +  GAMA(uj+1  +  uj_1 ) 

+  BETA(u.+1  j+1  -  u.+1  +  Ui  j_-,  -  ui_1  J+1)  +  i  UV  |  (4uJvl  -  uJV2) 

+  | VV |  (4uJvl  -  uJV2)]]}/{!  +  At[3(|UC|  +  j VC |  +  |UV|  +  |VV|) 

+  — -  (ALFA  +  GAMA)]}  (92) 

where  *  denotes  the  value  of  the  unknown  calculated  directly  from  the 
difference  equation. 

Several  observations  are  made  concerning  the  difference  equation  for 

u..  given  by  Equation  92.  The  unknowns  v..  and  p..  do  not  appear 
1  J  •  J  1  J 

explicitly.  In  fact  p . .  does  not  appear  at  all.  The  y  component  of 

*  J 

velocity  v..  does  appear  implicitly  in  the  quasi-1 inear  coefficients 

1  J 

UC  and  VC  as  does  u...  The  previous  value  for  u..  occurs  explicitly 

in  the  term  containing  the  velocity  divergence  D.  The  finite  difference 

linearization  assumption  requires  that  the  latest  prior  values  for  u.. 

■  J 

and  v^  are  used  in  these  instances. 
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In  an  analogous  manner,  the  y  component,  transformed  momentum  equation, 
(Equation  49)  is  written  for  .  and  becomes 

♦  att*ETA<pi+,  -  p,.,)  -  xxi(pj+1  -  Pj_,) 

+  |UC|  (4vJcl  -  vIC2)  +  I  VC  I  (4vjci  -  vJC2)  -  vD 
♦  pL  [ALFA(.1t,  ♦S»H»<vjtl  + 

+  BETA(v.+i  j+1  -  v.+]  J_1  +  v.^  -  v._i  j+1)  +  |UV|  (4vIvl  -  Vjy2) 

+  | VV|  (4vjy1  -  VJV2)]]}/{1  +  At [3 ( j UC |  +  |VC|  +  |UV|  +  |VV|) 

+  — -  (ALFA  +  GAMA]}  (93) 

Ret 

Again,  p..  does  not  appear  either  explicitly  or  implicitly.  The  unknown 

*  J 

u..  appears  implicitly  in  the  coefficients  (JC  and  VC.  Also,  v.  .  occurs 

'  J  1  J 

implicitly  in  UC  and  VC  and  explicitly  in  the  velocity  divergence  term. 

The  most  current  previous  values  for  u^  and  v^  are  used  in  these  lower 

order  coefficients. 

The  static  pressure  throughout  the  flow  field  is  calculated  from  the 
Poisson  pressure  equation  (Equation  50)  derived  in  Section  II.  Second 
order  accurate  central  differences  are  used  to  approximate  the  second 
derivatives  of  pressure  in  the  transformed  Laplacian  term.  The  first 
derivatives  of  pressure  in  the  Laplacian  are  approximated  by  second  order 
accurate  upwind  differences.  The  direction  is  based  on  the  sign  of  the 
coefficients  o  and  t  which  is  the  same  method  used  for  the  momentum  equa¬ 
tions.  In  this  equation,  a  nonlinear  source  term  occurs  which  is  composed 
of  transformation  and  velocity  derivatives.  The  transformation  derivatives 
are  approximated  with  second  order  accurate  central  differences.  Second 
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order  accurate  upwind  differences  are  used  for  the  velocity  derivatives. 
The  direction  of  the  differences  at  each  location  is  determined  by  the 
direction  of  the  local  £  and  n  velocity  components  related  to  UC  and  VC 
previously  defined.  The  time  derivative  of  the  velocity  divergence  is 
approximated  with  a  first  order  accurate  backward  difference.  The  require 
ment  that  the  velocity  divergence  at  the  nth  time  step  be  zero  is  incor¬ 
porated  into  the  difference  equation  by  setting  Dn  to  zero.  The  numerical 
non-zero  quantity  Dn  is  retained  as  a  correction  factor.  In  this  way, 
the  static  pressure  is  found  which  tends  to  satisfy  mass  conservation  at 
every  location  for  each  time  step. 

The  following  definitions  are  given  for  the  finite  difference  repre¬ 
sentation  of  the  source  terms  as  described  above: 

UX  .  (ynu5  -  ysun)/J 

uv  '  ,xe\  -  Vt>'J 

vx  -  <Vc  -  Vn)/J 

VY  =  (x_v  -  x  vr)/J 

£  n  n  £ 

If  these  definitions,  the  previous  transformation  coefficient  definitions, 
and  index  conventions  are  used,  the  transformed  pressure  equation  (Equa¬ 
tion  50)  is  rearranged  for  the  field  static  pressure  p? .  given  by 

J 

n  1 

pij  =  {*  ^At"  +  (UX)2  +  2(VX)  (UY)  +  (VY)2  +  ALFA  (pi+l  +  pi-l } 

+  GAMA(pj+1  +  Pj_-, )  +  BETA(pi+1  j+1  -  Pi+1  j_1  +  Pi_1  ^ 

-  Pi.!  j+1)  +  |UV|  (4pIy1  -  Piy2)  +  |W|  (4pjyl  -  Pjv2)>/(3(  |UV| 

+  | VV| )  +  2  (ALFA  +  GAMA)}  (94) 

In  this  equation  p. .  does  not  appear  implicitly.  The  velocity  components 
u^.  and  v^.  are  only  present  implicitly  in  the  nonlinear  source  terms. 
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The  three  difference  equations  (Equations  92  through  94)  form  a  set 
of  equations  for  the  flow  variables  u,  v,  and  p  at  location  (i,  j)  in  terms 
of  the  transformation  derivatives  and  values  of  the  variables  at  neighbor¬ 
ing  points.  The  quasi-linearization  and  differencing  methods  have  decoupled 
the  equations  with  respect  to  the  higher  order  terms  for  the  unknowns  at 
location  (i,  j). 


The  large  system  of  finite  difference  equations  which  must  be  solved 
for  the  n*"*1  time  step  is  solved  using  SOR  iteration  given  by 


»!**'>  «(;) 

ij  uv  ij  '  uv'  ij 

v^+1> ♦  0 .«  )  ,<;> 


1  j 


uv  IJ 


UV'  IJ 


P^+1)  =  (0  p*.  +  (1  -  «  )  P^} 

IJ  p  IJ  p  IJ 


(95a) 

(95b) 

(95c) 


where  s  is  the  iterate  number,  *  designates  the  current  value  obtained 
from  the  difference  equation,  and  a>uv  and  cup  are  acceleration  parameters. 

The  acceleration  parameters  are  always  unity  for  the  first  iteration. 

The  computation  of  the  wall  static  pressure  is  carried  out  separately 
as  discussed  in  Section  II.  The  pressure  iteration  technique  of  Chorin 
(Reference  68),  given  by  Equation  40,  avoids  the  difficulty  of  formulating 
one  sided  differences  at  the  body  for  the  Laplacian  of  pressure.  Also, 
at  the  body  surface,  conservation  of  mass  is  directly  imposed.  Second 
order  accurate  upwind  differences  are  used  for  the  transformation  and 
velocity  first  derivatives  in  the  velocity  divergence  of  Equation  51. 

The  acceleration  parameter  $,  which  is  related  to  SOR  iteration  of  the 
Poisson  pressure  equation,  is  given  by  Hodge  (Reference  53)  as 

4>  =  2Wpb/Ut(ALFA  +  GAMA)}  (96) 

where  u>pb  is  an  acceleration  parameter  for  the  corresponding  SOR  iteration 
of  the  Poisson  pressure  equation.  For  consistency,  the  iteration  equation 
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is  repeated  as 


(s) 

pb(s+l)  =  pb  -*D 


(97) 


where  s  represents  the  iterate  number  and  the  index  notation  is  kept  as 
described  above. 


The  trailing  edge  pressure  on  an  airfoil  can  be  difficult  to  obtain 
using  the  iteration  procedure  if  a  sufficiently  fine  grid  is  not  present. 
Hodge  (Reference  53)  experienced  such  difficulties  in  his  laminar  flow  work. 
The  trailing  edge  pressure  can  then  be  calculated  using  an  average  of  values 
obtained  by  linearly  extrapolating  the  pressures  calculated  at  the  closest 
points  on  the  upper  and  lower  surfaces.  For  a  two  point  linear  extra¬ 
polation  in  the  x  direction,  the  pressure  at  the  trailing  edge  becomes 
(x-,  -  x3)  p2  -  (Xl  -  x2)  p3 


2(x2  -  x3) 


(97a) 


XIMAX-2^  PIMAX-1  ‘  {xl  "xTMflv_i)  Pi 


IMAX-1'  pIMAX-2 


2  ( x 


IMAX-1 


XIMAX-2) 


where  the  x's  are  the  x  values  at  the  selected  surface  points  and  the 
surface  j  =  1  value  is  understood. 


The  iteration  procedure  follows  a  prescribed  order  at  each  point 
location  (i,  j).  Only  the  body  surface  pressure  is  computed  on  body  points. 
At  field  points,  the  u  velocity  component  is  computed  first,  then  the  v 
velocity  component,  and  finally  the  field  pressure  p.  The  points  are 
ordered  by  varying  the  n(j)  index  from  j  =  1  thru  j  =  0MAX-1  on  successive 
5  lines  in  the  transformed  plane.  The  5  lines  i  =  1  thru  i  =  IMAX-1  are 
traversed.  Values  for  the  variables  at  i  =  1  are  set  equal  to  the  values 
at  i  =  IMAX  to  ensure  continuity  across  the  transformation  branch  cut. 

The  i  index  is  adjusted  across  the  branch  cut  in  a  manner  identical  with 
that  described  for  obtaining  the  transformation  from  Equations  85  and  86. 
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3.  BOUNDARY  AND  INITIAL  CONDITIONS 

Boundary  and  initial  conditions  are  formulated  for  use  in  the  computa¬ 
tional  plane  to  completely  define  the  numerical  problem.  Boundary  condi¬ 
tions  on  the  body  surface  are  conveniently  specified  with  the  bodyfitted 
coordinate  system.  The  first  n  contour  is  constructed  to  be  the  body 
surface.  The  analytical  no  slip  boundary  conditions  for  velocity  given 
by  Equation  52  can  be  directly  applied  to  the  body  surface  in  the 
computational  plane  and  are  given  by 


ub  =  0 


vb  =  0 


where  the  subscript  b  denotes  the  body  surface  n  contour.  The  surface 
pressure  is  calculated  as  discussed  in  Section  IV. 2. 

The  5  lines  i  =  1  and  i  =  IMAX  in  the  computational  plane  define  the 
branch  cut  in  the  physical  plane.  No  boundary  conditions  are  required  or 
allowed  on  these  contours.  The  matching  of  solutions  on  these  lines,  which 
was  described  previously,  provides  a  continuous  solution  across  the  cut. 

The  outer  boundary  in  the  physical  plane  frequently  represents  con¬ 
ditions  at  infinity  (freestream  conditions)  given  by  Equations  53  and  54. 

In  the  body-fitted  coordinate  system,  the  outer  boundary  becomes  the  largest 
n  line  in  the  computational  plane.  The  physical  outer  boundary  is  chosen 
to  be  a  circle  with  a  radius  of  10  airfoil  chords.  The  distribution  of 
points  on  the  outer  boundary  is  readily  specified  on  a  circle;  also,  the 
initial  input  solution  for  obtaining  the  transformation  as  discussed  in 
Section  IV. 1  is  simple  to  specify.  Ghia  and  Hodge  (Reference  82)  applied 
a  numerical  inviscid  analysis  to  a  Joukowski  airfoil  at  angle  of  attack 
using  freestream  boundary  conditions  for  different  outer  boundaries. 

They  found  that  for  a  10°  angle  of  attack  the  lift  coefficient  differed 
by  less  than  one  percent  and  the  maximum  suction  pressure  coefficient  by 
about  two  percent  from  the  analytical  result.  Navier-Stokes  solutions 
using  a  circular  outer  boundary  with  a  smaller  radius  are  compared  with 
the  10  chord  radius  result  to  determine  effects  of  boundary  placement. 
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With  the  outer  boundary  elected  to  approximate  the  far  field  free- 
stream,  the  boundary  conditions  for  velocity  and  pressure  are  specified. 

As  pointed  out  by  Roache  (Reference  83),  caution  must  be  exercised  in 
applying  the  analytical  conditions  (Equations  53  and  54)  which  are  strictly 
valid  only  in  the  limit  of  large  distances  from  the  body.  If  these  condi¬ 
tions  are  used  on  a  boundary  at  a  finite  distance  away,  the  result  may 
predict  no  drag  since  a  wake  cannot  exist  at  the  outer  boundary.  This 
problem  is  avoided  by  applying  upstream  and  downstream  boundary  conditions 
on  the  outer  boundary.  The  upstream  boundary  along  which  conditions  are 
fixed  is  defined  by  the  semicircle  in  the  half  plane  x  ^  0  where  the  x 
axis  lies  along  the  airfoil  chord  with  the  origin  at  the  midchord.  The 
downstream  boundary  for  which  some  variables  are  permitted  to  be  free  is 
the  semicircle  in  the  half  plane  x  >  0. 

The  upstream  boundary  conditions  for  the  incoming  undisturbed  flow 
become 


ui  JMAX  C0S  a  vi  JMAX  "  sin  a  Pi  JMAX  =  0  (99) 

where  i  ranges  over  all  ?  contours  which  terminate  on  the  upstream  boundary. 

The  downstream  boundary  conditions  must  allow  a  wake  downstream  of  a 
body  to  pass  through  the  boundary.  The  no  change  boundary  condition  has 
been  widely  used  (References  53,83)  and  meets  this  requirement.  One 
approach  requires  no  change  of  the  velocity  components  in  the  freestream 
direction.  Using  the  expressions  for  the  directional  derivative  and 
gradient  in  the  transformed  plane,  this  boundary  condition  can  be  written 
for  the  u  component  of  velocity  as 

u  =  u,.  (x  sina  -  y  cosa)/  (xr  sina  -  y„  cosa)  (100a) 
n  5  n  n  5  £, 

where  a  is  the  geometric  angle  of  attack.  Now  consider  a  similar  down¬ 
stream  boundary  condition  where  the  no  change  condition  is  imposed  on  the 
velocity  components  along  the  downstream  £  contours.  In  a  similar  manner, 
this  condition  in  the  transformed  plane  becomes 

u  =  0  v  =  0  (100b) 

n  n 
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This  boundary  condition  is  identical  with  the  freestream  direction  boundary 
condition  of  Equation  100a  when  the  ?  contours  are  in  the  freestream  direction. 
In  this  case,  tana  =  y  /x  and  Equation  100a  becomes  identical  with  Equation 
100b.  This  case  does  occur  in  the  downstream  wake  region  for  the  trans¬ 
formation  obtained  in  Section  I\f.l.  Small  changes  in  the  velocity  compo¬ 
nents  are  found  in  the  far  wake  region  of  bodies.  The  downstream  far  field 
inviscid  flow  field  also  has  negligible  changes  as  seen  in  the  analysis 
of  Appendix  F.  Thus,  the  no  change  conditions  of  Equation  100b  approximate 
the  small  far  field  downstream  variations  in  velocity  while  allowing  the 
momentum  defect  in  a  wake  to  exit  at  the  outer  boundary. 

The  downstream  outer  boundary  condition  for  pressure  must  be  considered 
along  with  the  velocity  conditions.  Hodge  (Reference  53)  investigated  the 
use  of  a  similar  no  change  condition  for  pressure,  namely,  =  0  coupled 
with  the  velocity  conditions  of  Equation  100b.  His  laminar  solution 
developed  pressure  oscillations.  The  more  restrictive  condition  of 
specifying  the  freestream  pressure  in  the  far  wake  was  successfully 
used  by  Hodge  ^Reference  53).  The  specification  of  pressure  recovery 
with  a  freely  developing  velocity  defect  in  the  far  wake  gives  a  set  of 
numerical  boundary  conditions  which  can  physcially  represent  the  far  wake 
at  a  finite  far  field  downstream  boundary.  The  downstream  pressure 
boundary  condition  becomes 


pi  JMAX  "  0  (101) 

The  suitable  description  of  the  physics  at  the  downstream  computa¬ 
tional  boundary  must  be  written  in  finite  difference  form.  The  no  change 
condition  approximated  with  a  second  order  accurate  central  difference 
applied  at  the  midpoint  in  the  computational  plane  between  the  boundary 
point  and  first  interior  point  on  a  downstream  £  contour  (i  index)  gives 

UJMAX  =  UJMAX-1  V JMAX  =  VJMAX-1  *102^ 


for  the  downstream  boundary  conditions  of  the  velocity  components.  An 
alternate  formulation  is  obtained  by  approximating  the  no  change  condition 
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with  a  second  order  accurate  upwind  difference  applied  at  the  boundary 
point.  Then,  the  downstream  boundary  conditions  for  the  velocity  components 
become 


UJMAX  =  ^4uJMAX-1  "  UJMAX-2^/3 
VJMAX  =  (4vJMAX-1  "  vJMAX-2^/3 


(103) 


These  boundary  conditions  also  follow  from  a  second  degree  polynominal  fit 
with  a  zero  derivative  imposed  at  the  downstream  boundary. 


The  time  dependent  finite  difference  method  given  previously  requires 
a  set  of  initial  values  for  the  velocities  and  pressures  at  all  the  (i,  j) 
locations.  This  requirement  parallels  the  result  discussed  for  the  dif¬ 
ferential  equations  in  Section  III. 2.  Two  methods  are  used  to  provide 
initial  conditions.  In  the  first  method,  an  inviscid  solution  with  cir¬ 
culation  for  the  flow  field  at  the  given  angle  of  attack  is  obtained  using 
SOR  iteration  as  described  in  Appendix  G.  This  solution  then  provides 
initial  values  for  the  Navier-Stokes  solution.  The  second  technique  uses 
the  Navier-Stokes  solution  for  one  angle  of  attack  and  rotates  the  velocities 
by  the  change  in  angle  of  attack.  The  set  of  initial  values  for  velocity 
and  pressure  is  then  used  to  start  a  Navier-Stokes  solution  at  the  new 
angle  of  attack. 


4.  TURBULENCE  MODEL 

The  numerical  procedure  which  implements  the  eddy  viscosity  turbulence 
model  described  in  Section  III  is  discussed.  The  two-layer  model  near  the 
airfoil  surface  is  presented  followed  by  the  far  wake  model.  The  procedure 
used  in  the  near  wake  is  then  described.  The  eddy  viscosity  distribution 
throughout  the  flow  field  is  calculated  at  the  beginning  of  each  time  step. 
SOR  iteration  is  then  used  to  determine  the  velocities  and  pressures  at  the 
new  time  level .  The  eddy  viscosity  distribution  is  not  changed  during  the 
iteration  process. 
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The  inner  viscosity  computation  from  Equation  60  requires  values  of 
the  tangential  and  normal  velocity  derivatives  and  the  distance  from  the 
surface.  These  velocity  derivatives  are  given  by  Equations  70  and  71  in 
the  transformed  plane.  They  are  evaluated  using  second  order  accurate 
central  differences  for  both  the  transformation  derivatives  and  the  u  and  v 
velocity  component  derivatives  with  respect  to  £  and  n.  The  inner  eddy 
viscosity  is  identically  zero  at  the  surface  (n=l)  and  need  not  be 
calculated  separately.  The  corresponding  velocity  expression  which  is 
found  in  the  exponent  of  the  damping  factor  (Equation  62)  must  be  evaluated 
at  the  airfoil  surface  for  each  £  contour  in  regions  of  turbulence.  The 
derivatives  with  respect  to  n  are  approximated  by  second  order  accurate 
forward  differences  found  in  Appendix  B,  and  the  derivatives  with  respect 
to  £  are  evaluated  using  second  order  accurate  central  differences.  The 
distance  from  the  surface  to  the  point  (i,  j)  is  measured  along  the 
constant  £  line  represented  by  index  i.  The  freestream  chord  Reynolds 
number  appears  throughout  the  model  as  a  result  of  writing  the  equations 
in  nondimensional  form  and  is  an  input  parameter. 

The  adverse  pressure  gradient  turbulence  model  modification  in  the 
trailing  edge  region  given  by  Equation  65  is  implemented  next.  The  origin 
of  the  distance  s,  measured  along  the  airfoil  surface  from  the  trailing  edge 
separation  point,  is  determined  by  examining  the  u  velocity  component  at 
points  on  the  n  contour  next  to  the  surface  (j=2).  The  search  begins  on 
the  £  line  corresponding  to  the  95  percent  chord  location  and  advances 
toward  the  airfoil  leading  edge  along  the  second  n  line.  The  search 
terminates  when  a  positive  value  for  u.£  is  detected  which  indicates 
attached  flow.  The  relaxation  factor  '  is  then  evaluated  for  each 
applicable  £  line  where  trailing  ”do  ^paration  is  detected.  The 
inner  eddy  viscosity  is  multip'  .■  oy  .  ■  relaxation  factor  which 
adjusts  the  adverse  pressure  gradient  inner  eddy  viscosity  constant  k-j . 
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After  the  inner  eddy  viscosity  is  computed,  the  limiting  technique 
evaluates  the  computed  value.  The  inner  eddy  viscosity  is  prevented  from 
decreasing  in  the  outward  normal  and  downstream  tangential  flow  directions 
by  comparing  the  value  with  previously  computed  values  as  given  by  the 
following  sequence; 


e .  .  >  e .  .  ,  and 
1J  -  l  J-l 


1J  -  1-1  J 


The  larger  value  is  selected  during  each  comparison.  The  initial  inner 
eddy  viscosity  limiting  distribution  is  calculated  in  the  attached  boundary- 
layer  region  with  a  favorable  pressure  gradient  near  the  airfoil  leading 
edge.  The  first  or  outward  direction  limiting  condition  is  imposed  here 
as  well . 


The  limit  value  of  the  inner  eddy  viscosity  is  then  compared  with  the 
calculated  outer  layer  eddy  viscosity.  Whenever  the  inner  eddy  viscosity 
first  exceeds  the  outer  value,  the  outer  eddy  viscosity  is  then  used  for  the 
remaining  locations  outward  from  the  surface.  This  switch  ensures  the 
continuity  of  the  eddy  viscosity  distribution. 

The  computation  of  the  outer  eddy  viscosity  given  by  Equation  61 
requires  values  for  the  local  boundary-! ay er  edge  velocity  and  the 
displacement  thickness  defined  by  Equation  64.  These  quantities  are 
calculated  for  each  ?  contour  which  crosses  the  turbulent  boundary-layer. 

For  a  c  line,  the  local  edge  velocity  ufi  is  defined  to  be  the  maximum 
tangential  velocity  measured  relative  to  the  local  surface  tangent  line. 

The  tangential  velocity  is  given  by  Equation  68  and  is  evaluated  as 
described  previously  for  the  inner  eddy  viscosity.  The  boundary-layer 
thickness  6  is  then  the  distance  along  the  5  line  from  the  surface  to 
the  first  point  where  the  tangential  velocity  is  not  less  than  0.99  ug. 

The  displacement  thickness  is  calculated  using  trapezoidal  integration 
where  the  limits  of  integration  are  determined  by  the  computed  boundary- 
layer  thickness. 
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The  resulting  eddy  viscosity  distribution  is  modified  next  to 
simulate  the  transition  from  laminar  to  fully  turbulent  flow.  The 
transition  factor  r,  given  by  Equation  67,  is  calculated  as  follows. 

The  starting  location  for  transition  on  the  airfoil  surface  must  be 
specified.  This  location  is  specified  by  designating  the  first  grid 
point  on  the  surface  downstream  as  well  as  the  distance  to  that  first 
grid  point  as.  The  transition  distance  s  from  the  starting  location 
to  a  given  £  line  is  calculated  by  adding  as  to  the  distance  between 
the  successive  surface  grid  points  with  known  locations.  The  relaxation 
factor  A  is  determined  from  the  transition  factor  Equation  67  and  the 
definition  of  A.  If  r  =  3/4  is  specified  when  the  transition  distance 
equals  three  fourths  of  the  total  transition  length,  then  1/A  is  given 
by 


1/A  =  2.4456/xl  (104) 

where  xL  is  the  total  assumed  transition  length.  The  eddy  viscosity 
distribution  computed  previously  on  a  £  line  is  multiplied  by  the 
corresponding  constant  r  to  obtain  the  revised  distribution  of  turbulent 
eddy  viscosity  which  includes  the  transition  region. 

The  final  modification  made  to  the  eddy  viscosity  distribution  across 
the  upper  surface  turbulent  region  involves  the  decrease  of  turbulence  in 
the  direction  from  the  boundary-layer  toward  the  far  field.  The  inter- 
mittency  factor  given  by  Equation  63  models  this  behavior.  The  previously 
calculated  boundary-layer  thickness  5  and  the  distance  from  the  airfoil 
surface  to  the  (i,  j)  location  on  a  5  line  are  used  to  compute  the 
intermittency  y  for  that  location.  The  previous  values  of  the  outer 
eddy  viscosity  are  multiplied  by  the  corresponding  value  of  the  inter¬ 
mittency  factor  to  obtain  the  final  eddy  viscosity  distribution  for  the 
turbulent  region  above  the  airfoil  surface  at  the  nth  time  step. 
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The  procedure  for  computing  the  eddy  viscosity  distribution  on  a  £ 
line  across  the  turbulent  flow  zone  above  the  airfoil  is  repeated  down¬ 
stream  for  each  c  line  through  the  95  percent  chord  point.  The  £  contours 
near  the  trailing  edge  deviate  from  the  local  normal  direction  in  excess 
of  10°.  For  locations  in  this  small  region,  the  eddy  viscosity  is  found 
by  using  the  last  calculated  value  upstream  in  a  direction  tangent  to  the 
airfoil  surface. 

The  calculation  of  the  far  wake  eddy  viscosity  given  by  Equation  73 
is  similar  to  the  computation  of  the  outer  layer  eddy  viscosity.  In  the 
wake  region  the  n  contours  are  approximately  orthogonal  to  the  flow 
direction  in  the  wake.  The  n  contours  are  thus  conveniently  defined  paths 
across  the  wake  for  use  in  computing  wake  characteristics. 

On  an  n  contour,  the  component  of  velocity  in  the  freestream  direction 
is  computed  for  the  grid  points  spanning  the  wake.  The  minimum  of  these 
directional  velocity  components  in  the  wake  is  determined.  The  maximum 
values  of  the  velocity  components  which  occur  on  either  side  of  the  minimum 
value  are  found  and  designated  the  upper  and  lower  edge  velocities.  Then 
the  upper  and  lower  edges  of  the  wake  are  defined  as  the  points  where  the 
local  velocity  component  in  the  freestream  direction  first  reaches  99  per¬ 
cent  of  the  upper  and  lower  edge  velocities,  respectively.  The  distances 
from  the  point  of  minimum  velocity  to  the  upper  and  lower  wake  edges  are 
defined  as  the  upper  and  lower  wake  thicknesses,  respectively.  These 
distances  determine  the  limits  of  integration  when  calculating  the  upper 
and  lower  wake  displacement  thicknesses  from  Equation  64.  The  origin 
of  the  variable  of  integration  is  at  the  grid  point  where  the  minimum 
velocity  occurs.  The  upper  and  lower  wake  displacement  thicknesses  are 

numerically  evaluated  with  trapezoidal  integration.  The  half  wake  dis- 

★ 

placement  thickness  6  is  the  average  of  the  upper  and  lower  displacement 

w 

thicknesses.  The  far  wake  eddy  viscosity  is  calculated  from  Equation  73 

★ 

using  the  thickness  6  and  an  average  of  the  upper  and  lower  edge 
velocities . 
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The  intermittency  factor  for  each  grid  point  on  an  n  contour  is 
determined  using  the  average  value  of  the  upper  and  lower  wake  thicknesses 
for  6  in  Equation  63.  The  local  distance  to  each  point  is  calculated 
from  the  location  of  minimum  velocity  along  the  n  contour. 

The  near  wake  model  is  numerically  implemented  by  considering  the 
upper  turbulent  and  lower  laminar  boundary-layers  separately.  The  eddy 
viscosity  distribution,  which  is  calculated  near  the  trailing  edge  of  the 
airfoil  at  the  95  percent  chord  location,  is  extended  into  the  near  wake 
along  a  direction  parallel  to  the  airfoil  surface  defined  by  the  95  percent 
chord  and  99  percent  chord  grid  point  locations.  The  chord  line  is  the 
line  of  extension  for  the  boundary-layers  in  the  near  wake.  The  £  lines 
i  =  1  and  i  =  IMAX-1  define  the  interaction  zone  in  the  near  wake.  The 
noninteracting  distance  WL-j  in  Equation  74  is  based  on  the  length  of  the 
lower  laminar  boundary-layer  thickness  computed  at  the  95  percent  chord 
location.  This  boundary-layer  thickenss  is  computed  by  the  same  method 
used  previously  for  the  upper  surface  boundary-layer.  The  distance  WLZ 
along  the  chord  line  from  the  trailing  edge  to  the  location  where  the 
far  wake  eddy  viscosity  model  begins  is  also  calculated  as  a  multiple 
of  the  laminar  trailing  edge  boundary-layer  thickness.  The  eddy  viscosity 
in  the  interaction  zone  is  computed  with  the  two  calculated  distances  and 
the  far  wake  eddy  viscosity  using  Equation  74. 

5.  FORCE  COEFFICIENTS 

The  forces  and  moments  exerted  on  the  airfoil  body  by  the  flow  field 
are  evaluated  at  the  body  surface  using  Equations  79,  80,  and  82.  The 
integrals  are  evaluated  using  trapezoidal  integration.  The  product  terms 
in  the  integrands  of  the  form  (fg)  are  integrated  by  using  products  of 
averages,  (f..  +  f..+^)  (g^  +  +-j  )/4 .  The  £  derivatives  are  evaluated 

using  second  order  accurate  central  differences,  and  the  n  derivatives 
are  formulated  with  second  order  accurate  upwind  forward  differences. 

The  lift  and  drag  coefficients  are  calculated  by  adding  the  components 
of  x  and  y  direction  force  coefficients  in  the  normal  freestream  (90 
degrees  counterclockwise)  and  freestream  directions,  respectively,  using 
Equations  D-13  and  D-14. 
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The  alternate  technique  based  on  the  control  volume  analysis  in 
Appendix  D  for  computing  the  forces  on  the  airfoil  body  is  numerically 
implemented  in  a  similar  manner.  The  force  coefficients  in  the  x  and  y 
directions  for  a  two-dimensional  control  volume  defined  by  an  n  line  are 
given  by  Equations  83  and  84,  respectively.  The  line  integrals  are  again 
computed  with  trapezoidal  integration  using  products  of  averages.  The  E 
derivatives  are  evaluated  with  second  order  accurate  central  differences. 

The  n  derivatives  are  evaluated  with  second  order  accurate  central  dif¬ 
ferences  on  the  interior  n  lines.  Second  order  accurate  upwind  differences 
are  used  to  approximate  the  n  derivatives  on  the  body  and  outer  boundary 
n  lines.  The  area  integral  is  also  evaluated  using  trapezoidal  integration. 
The  Jacobian  for  an  area  element  bounded  by  the  i  and  i+1  E  lines  and  j-1 
and  j  n  lines  is  calculated  by  taking  the  average  of  the  Jacobian  values 
previously  computed  for  the  viscous  term  in  the  line  integrals  at  the  cor¬ 
responding  E  segments  on  the  j-1  and  j  n  lines.  The  area  integrals  are 
computed  for  both  n  and  n-1  time  intervals  using  the  previously  computed 
flow  field  values.  The  time  derivative  is  then  evaluated  using  a  first 
order  accurate  backward  difference. 
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SECTION  V 

DISCUSSION  OF  RESULTS 

Numerical  solutions  for  two-dimensional  incompressible  turbulent 
viscous  flow  over  a  NACA  0012  airfoil  section  are  obtained  using  the 
implicit  finite  difference  method  previously  described.  The  NACA  0012 
airfoil  section  is  chosen  because  of  its  widespread  use  as  a  test  case 
in  both  experimental  and  computationa1  work  (References  65,66).  The 
solutions  are  for  angles  of  attack  of  5°,  7.5°,  9.5°,  and  11.5°  at  a 
chord  Reynolds  number  of  170,000.  This  Reynolds  number  is  selected  for 
comparison  with  both  previous  NACA  data  (Reference  56)  and  Air  Force 
F.  J.  Seiler  Research  Laboratory  data  (Reference  84). 

The  numerically  generated  body-fitted  coordinate  system  which  was 
used  for  each  solution  is  presented.  The  successive-over-relaxation 
(SOR)  iteration  numerical  parameters  and  convergence  criteria  along  with 
the  steady  state  convergence  criteria  are  then  discussed.  The  velocity 
and  pressure  fields  and  streamline  contours  obtained  from  the  numerical 
solutions  are  examined.  The  calculated  laminar  separation  bubble  charac¬ 
teristics  are  compared  with  empirical  results.  The  position  of  transition 
to  turbulence  and  the  transition  length  relative  to  the  separation  bubble 
location  are  discussed.  The  computed  airfoil  surface  mean  pressure  distri¬ 
butions  are  compared  with  both  an  inviscid  solution  and  experimental  data. 
The  numerically  predicted  lift  and  drag  forces  exerted  by  the  flow  field  on 
the  airfoil  are  presented  along  with  available  experimental  measurements. 
The  effects  that  the  placement  and  type  of  imposed  far  field  boundary 
conditions  have  on  the  numerical  solution  are  discussed.  The  sensitivity 
of  the  numerical  solution  to  grid  size  is  examined. 

The  curvilinear  body-fitted  grid  shown  in  Figure  3  is  numerically 
generated  using  boundary-layer  parameters  for  the  attraction  of  n  con¬ 
tours  towards  the  airfoil  surface  and  linear  interpolation  for  wake 
resolution  given  by  Appendix  C.  The  Blasius  series  flat  plate  chord 
Reynolds  number  of  200,000,  a  nondimensional ized  boundary-layer  edge 
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velocity  of  1.2,  and  seven  grid  points  specified  in  the  boundary-layer 
are  used  in  the  n  contour  attraction  technique  (Reference  7).  The 
distribution  of  airfoil  surface  points  is  accomplished  with  Equation  26 
where  A1  =  0.6,  A2  =  0.4,  and  A3  =  1.0.  Seventy-one  6  or  surface  points 
together  with  44  n  contours  are  used.  The  transformation  difference  equa¬ 
tions  (Equations  85  and  86)  were  then  solved  using  SOR  iteration.  The 
optimum  local  acceleration  parameters  varied  between  0.8  and  1.5  throughout 
the  flow  field  with  essentially  no  change  in  the  local  values  after  25 
iterations.  After  350  iterations,  the  maximum  errors  in  the  solution 
for  both  x  and  y  were  less  than  |10~^|  and  occurred  near  the  outer  boundary 
above  the  airfoil.  Additional  6  lines  were  added  in  the  airfoil  wake 
region  using  the  interpolation  method  described  in  Appendix  C.  Eight 
grid  points  (four  at  a  time)  were  inserted  on  the  rounded  trailing  edge. 
Table  C-l  gives  the  6  and  n  indexing  for  the  final  79x44  grid  shown  in 
Figure  3  and  the  equivalent  indices  for  the  71x44  original  grid. 

The  airfoil  is  thus  defined  by  a  total  of  79  points.  Eleven  grid 
points  define  the  first  5%  chord  nose  region  and  17  points  are  used  to 
define  the  rounded  trailing  edge.  The  maximum  distance  between  successive 
surface  grid  points  is  5%  of  the  chord.  The  6  lines  intersect  the  surface 
n  contour  within  10°  of  the  local  normal  direction  except  in  the  last  5% 
chord  region  near  the  trailing  edge. 

The  SOR  iteration  parameters  and  convergence  criteria  which  are 
required  by  the  numerical  implicit  procedure  are  similar  for  each  angle 
of  attack  solution.  The  acceleration  parameters  had  values  of  1.0  for  the 
velocity  components,  0.9  for  the  surface  pressure,  and  1.10  or  1.15  for 
field  pressures.  The  larger  value  of  was  used  at  9.5°  and  11.5°  angles 
of  attack.  The  convergence  criteria  for  the  SOR  iteration  required  that 
the  maximum  change  in  the  relative  magnitudes  of  u,  v,  and  p  for  all 
locations  be  less  than  1%,  This  criteria  was  relaxed  to  5%  for  the  9.5° 
and  11.5°  cases.  Three  iterations  per  time  step  were  the  maximum 
necessary  during  most  of  each  computed  case.  For  the  9.5°  case,  the 
approach  to  the  steady  state  solution  was  re-run  with  <o  ^  =  0.95  and  Wp  = 

1.10  for  two  characteristic  times.  An  additional  iteration  per  time  step 
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was  required  to  maintain  the  same  error  tolerance.  Computed  values  of 
lift  and  drag  were  virtually  identical. 

Although  implicit  SOR  methods  have  been  shown  by  linear  stability 
analyses  to  be  unconditionally  stable,  convergence  at  a  given  time  re¬ 
quired  a  small  time  step  At  <  0.001.  The  analysis  in  Appendix  E  based 
on  linear  matrix  theory  indicates  that  the  diagonal  dominance  of  the  set 
of  finite  difference  equations  given  by  Equations  92  and  93  has  a  definite 
dependence  on  1/At.  Diagonal  dominance  is  a  necessary  and  sufficient 
condition  for  convergence  (Reference  81)  of  constant  coefficient  linear 
systems  using  SOR  iteration.  Although  this  system  of  equations  has  lower 
order  nonlinearities  and  variable  coefficients,  the  order  of  magnitude 
study  of  Appendix  E  indicates  convergence  if  the  most  stringent  require¬ 
ment  At  <  0.0003  is  applied.  A  constant  time  step  equal  to  0.0005  was 

used  for  each  solution.  No  convergence  problems  were  encountered. 

The  finite  difference  method  described  in  Section  IV  is  a  time  de¬ 
pendent  technique.  The  solution  is  advanced  in  time  by  At  increments. 

The  numerical  solution  was  considered  to  have  reached  a  steady  state  when 
changes  in  the  computed  values  of  the  lift  and  drag  coefficients  were  less 
then  1%  over  one  characteristic  or  nondimensional  time  period.  Each  steady 
state  solution  required  three  to  six  characteristic  time  periods  and  from 
two  to  four  hours  of  CPU  time  on  a  CDC  CYBER  750  computer. 

The  numerical  solution  of  the  mean  flow  field  near  the  airfoil  sec¬ 
tion  at  each  angle  of  attack  is  presented  by  using  streamline  contours, 
velocity  field  vectors,  and  pressure  contours.  The  mean  streamline  con¬ 
tours  are  shown  in  Figures  4  through  7.  A  small  laminar  separation  bubble 
on  the  upper  surface  with  negligible  trailing  edge  separation  is  seen  in 
Figure  4  for  a  =  5°.  For  an  increased  angle  of  attack  equal  to  7.5°,  the 
separation  bubble  has  decreased  in  size  and  moved  forward  toward  the  leading 
edge  of  the  airfoil.  Figure  5  also  shows  the  turbulent  trailing  edge 
separation  region  which  has  also  moved  forward.  The  trailing  edge 
separation  region  has  grown  significantly  and  has  progressed  to  the 
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quarter  chord  position  when  the  angle  of  attack  has  reached  11.5°.  This 
behavior  is  typically  observed  during  trailing  edge  stall  (Reference  72). 

The  velocity  field  vector  plots  presented  in  Figures  8  through  11 
also  illustrate  the  phenomena  discussed  previously.  In  addition,  a  small 
scale  clockwise  rotational  motion  is  observed  in  the  rear  portion  of  the 
trailing  edge  separation  region.  A  circulatory  flow  pattern  within  the 
laminar  separation  bubble  is  seen  for  a  =  5°  in  Figure  8.  The  boundary- 
layer  near  the  leading  edge  for  each  angle  of  attack  contains  at  least 
five  grid  points  along  each  £  contour  while  approximately  10  points  resolve 
the  lower  surface  laminar  boundary-layer  near  the  trailing  edge.  This 
number  of  points  should  be  sufficient  to  resolve  the  primary  features. 

The  pressure  contours  near  the  airfoil  in  Figures  12  through  15  clearly 
show  expanding  regions  of  low  pressure  above  the  airfoil  and  high  pressure 
below  the  airfoil  as  the  angle  of  attack  increases.  The  large  pressure 
variations  in  the  upper  surface  suction  peaks  are  also  observed.  The 
pressure  fields  indicate  that  only  small  pressure  variations  occur  in 
the  wake  for  each  angle  of  attack.  The  smooth  and  continuous  numerical 
solution  for  the  pressure  field  is  in  sharp  contrast  to  the  large 
oscillatory  results  obtained  by  Hodge  (Reference  7)  for  laminar  flow. 


For  a  moderately  thick  airfoil  with  a  smooth  surface  in  a  flow 
with  a  low  freestream  turbulence  intensity,  the  laminar  boundary-layer 
near  the  leading  edge  can  separate  upon  encountering  a  strong  adverse 
pressure  gradient.  Subsequently,  shear  layer  instabilities  can  cause 
transition  to  turbulence.  The  increased  mixing  may  reattach  the  shear 
layer  and  thereby  form  a  separation  bubble.  The  formation  of  separation 
bubbles  and  their  relationship  witn  the  stalling  characteristics  of 
airfoils  at  various  Reynolds  numbers  have  been  investigated  by  Gault 
(Reference  85),  Gaster  (Reference  86),  Hoad,  et  al  (Reference  87),  and 
Arena  and  Mueller  (Reference  88).  The  laminar  separation  bubble  has  been 
observed  experimentally  by  Gregory  and  O'Reilly  (Reference  89)  for  the 
NACA  0012  airfoil  over  a  large  range  of  Reynolds  numbers.  These  experi¬ 
mental  and  empirical  results  will  be  used  to  compare  with  the  numerical 
solutions . 
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In  this  investigation,  the  beginning  of  turbulent  transition  on 
the  upper  surface  of  the  airfoil  and  the  transition  length  have  been  based 
on  the  closure  of  the  laminar  separating  shear  layer.  The  surface  mean 
pressure  coefficients  for  a  =  5°  and  7.5°  presented  in  Figures  16  and 
17  clearly  show  the  laminar  bubble  constant  pressure  region,  downstream 
of  the  suction  peak,  followed  by  a  rapid  recovery.  Transition  in  the 
numerical  solution  occurs  near  the  downstream  side  of  the  bubble  where 
the  steep  pressure  recovery  begins.  This  phenomenon  has  been  reported  by 
Wallis  (Reference  90}  and  Arena  and  Mueller  (Reference  88).  The  chord 
locations  for  the  start  of  transition  xt  and  full  turbulence  x^  relative 

to  the  separation  and  reattachment  points  of  the  bubble,  x$  and  xR 

respectively,  are  given  in  Table  1.  For  the  smaller  angles  of  attack, 
reattachment  occurred  prior  to  attaining  fully  turbulent  flow.  This 
behavior  has  been  experimentally  observed  by  Arena  and  Mueller  (Reference 
88)  in  their  low  Reynolds  number  flow  study. 

The  computed  values  of  several  parameters  at  separation  for  both 
a  =  5°  and  a  =  7.5°  cases  agree  with  empirical  results  and  are  presented 
in  Table  2.  The  pressure  gradient  parameter  M  =  -  B2  Re  dug/ds  in 
nondimensional  variables,  where  0  is  the  momentum  thickness  and  s  is 
arclength,  compares  favorably  with  the  laminar  separation  criteria  of 
Curie  and  Skan  (Reference  91)  given  by  M  0.09.  The  shape  parameter 
H12  =  <S*/,0  tal<es  on  va^ues  of  3.1  and  3.2  at  laminar  separation  compared 
with  an  average  value  of  3.5  reported  by  Curie  and  Skan.  The  Reynolds 
number  at  separation  based  on  the  displacement  thickness  (Re{*s)  has  a 

value  greater  than  500  and  predicts  a  short  bubble  using  the  Owen-Klanfer 
(Reference  92)  criteria  confirmed  by  Gault  (Reference  85)  and  Crabtree 
(Reference  93).  The  bubble  length  Lg  is  of  the  order  10^  s*$  and  decreases 

rapidly  with  increased  angle  of  attack  as  seen  in  Table  2.  This  relation¬ 
ship  also  indicates  a  short  bubble  (References  90,93)  whereas  long  bubbles 
have  lengths  of  the  order  10^  5*^  (References  85,86).  The  assumption 

that  laminar  separation  precedes  turbulent  transition  is  verified  with 
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TABLE  1 

COMPUTED  CHARACTERISTICS  OF  TURBULENCE  NEAR 
THE  AIRFOIL  LEADING  EDGE  (x  =  -0.5) 


a  (Deg) 

XS 

Xt 

XR 

XT 

5.0 

-  0.42 

-  0.24 

-  0.15 

-  0.09 

7.5 

-  0.47 

-  0.42 

-  0.38 

-  0.35 

9.5 

- 

-  0.48 

- 

-  0.46 

11.5 

_ 

-  0.49 

_ 

-  0.48 

TABLE  2 

COMPUTED  LAMINAR  SEPARATION  BUBBLE  CHARACTERISTICS 


a  (Deg) 

HS 

H12S 

Re6*S 

ReeS 

S* 

S 

lb 

AR 

5 

.124 

3.1 

680 

220 

.0025 

.27 

-.002 

7.5 

.091 

3.2 

540 

168 

.0013 

.09 

CM 

o 

o 

1 
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the  computed  Reynolds  number  at  separation  based  on  momentum  thickness 
Re0g.  Crabtree's  (Reference  94)  criteria  indicates  that  transition  has 

occurred  if  ReQ^  >  780  when  first  reaches  0.09.  Thus,  the  solutions 

correctly  predict  that  laminar  separation  occurs  before  transition  to 
turbulence.  The  turbulent  reattachment  criterion  for  the  bubble  given  by 
Roberts  (Reference  95)  is  A0  =  (©  /u  )  du  /ds  >  -  0.0059.  This  criterion 

is  also  satisfied  as  shown  in  Table  2. 

The  trends  discussed  previously  for  M^,  A^,  and  both  boundary-layer 

Reynolds  numbers  continue  at  the  higher  angles  of  attack.  However,  at 
angles  of  attack  near  stall  the  bubble  length  decreases  to  about  1%  chord 
(Reference  90)  which  is  the  order  of  the  £  contour  spacing  near  the  leading 
edge.  The  numerical  solutions  at  a  =  9.5°  and  a  =  11.5°  consequently  do 
not  resolve  the  bubble,  as  seen  in  Figures  18  and  19.  The  solutions 
become  very  sensitive  to  small  changes  in  the  turbulent  transition  length 
xL  and  starting  location  x^.  For  example,  at  a  =  9.5°  with  an  initial 

transition  location  of  xt  =  -0.48,  a  change  in  transition  length  from 
0.022  to  0.025  produced  shedding  of  bubbles,  large  pressure  oscillations, 
and  a  significant  30%  increase  in  average  drag  compared  with  the  steady 
state  solution.  The  experimental  pressure  data  were  measured  with  dia¬ 
phragm  transducers  capable  of  measuring  frequencies  well  beyond  5000  Hz. 

No  dominant  frequency  was  observed.  Hoad,  et  al  (Reference  87)  recently 
observed  a  small  scale  unsteadiness  in  a  short  region  near  the  surface 
of  a  NACA  0012  airfoil  at  angle  of  attack.  The  laser  velocimeter  mean 
velocity  histograms  for  a  certain  region  near  the  leading  edge  exhibited 
dominant  and  small  secondary  mean  velocities.  Hoad,  et  al  suggest  that 
this  behavior  may  indicate  an  unsteadiness  within  the  laminar  separation 
bubble.  However,  the  remaining  flow  field  was  steady. 

The  numerical  solutions  presented  for  each  angle  of  attack  are  the 
result  of  a  parametric  study  where  the  upper  surface  transition  length  and 
location  were  varied  to  close  the  bubble  and  maintain  steady  flow.  In 
each  case,  a  further  increase  in  the  transition  length  or  movement  down¬ 
stream  of  the  start  of  transition  resulted  in  a  nonphysical  oscillatory 
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motion  propagated  downstream  with  an  accompanied  large  increase  in  drag. 

The  good  agreement  between  the  calculated  and  empirical  separation  bubble 
characteristics  indicates  that  the  origin  of  turbulence  on  the  upper  sur¬ 
face  is  satisfactorily  modelled. 

Computed  boundary-layer  velocity  profiles  at  four  stations  on  the 
upper  surface  are  compared  with  two  sets  of  hot  wire  anemometry  data 
(Reference  84)  in  Figures  20  through  23.  The  profiles  were  measured 
and  computed  at  the  same  locations  on  constant  E,  lines  located  approxi¬ 
mately  at  18%,  50%,  69%,  and  92%  chord.  The  profiles  are  nondimensional ized 
by  the  locally  computed  boundary-layer  thickness  which  is  defined  as  the 
normal  distance  from  the  surface  to  where  the  tangential  velocity  first 
attains  99%  of  the  local  edge  velocity.  Significant  RMS  fluctuations 
for  each  angle  of  attack  were  detected  experimentally  at  the  inner  data 
point  locations  for  stations  two  through  four  which  indicate  the  presence 
of  a  turbulent  boundary-layer.  No  fluctuations  were  observed  in  the  lower 
surface  boundary-layer  at  the  90%  chord  location.  This  result  provides 
experimental  evidence  for  the  assumption  that  the  lower  surface  boundary- 
layer  remains  laminar  in  the  presence  of  the  favorable  pressure  gradient 
at  the  Reynolds  number  under  investigation.  The  calculated  boundary-layer 
becomes  thicker  downstream  compared  with  the  experimental  profiles.  This 
result  could  be  caused  by  grid  boundary-layer  resolution  errors  propagated 
downstream  or  by  deficiencies  in  the  turbulence  model.  The  hot  wire  data 
have  an  estimated  experimental  error  of  5%. 

The  boundary-layers  on  the  upper  and  lower  surfaces  merge  to  form  the 
wake  downstream  of  the  airfoil.  Wake  mean  velocity  profiles  at  chord 
locations  of  0.58,  0.79,  and  1.54  for  each  angle  of  attack  are  shown  in 
Figures  24  through  27.  The  trailing  edge  is  located  at  x  =  0.5.  The 
profiles  are  measured  along  constant  n  lines  with  the  origin  on  the 
extended  chordline.  The  thicker  numerical  upper  surface  boundary-layer 
propagates  downstream.  The  computed  large  variations  of  velocity  near 
the  lower  edge  of  the  wake,  which  come  from  the  lower  surface  boundary- 
layer,  are  in  good  agreement  with  the  hot  wire  anemometry  measurements. 
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The  measurements  were  obtained  at  grid  points  using  the  original  71x44 
grid  while  the  numerical  solutions  were  obtained  from  the  79x44  grid  system 
described  previously.  The  computed  values  of  the  wake  displacement  and 
momentum  thicknesses  approach  approximately  constant  values  within  15% 
of  each  other  at  a  location  of  one  chord  downstream  of  the  trailing  edge. 
They  were  computed  using  second  order  accurate  trapezoidal  integration 
and  are  given  in  Table  3. 


TABLE  3 

COMPUTED  FAR  WAKE  CHARACTERISTICS 


a  (Deg) 

5* 

w 

Qw 

5.0 

0.011 

0.010 

7.5 

0.019 

0.017 

9.5 

0.028 

0.026 

11  .5 

0.038 

0.033 

_  2 

The  Reynolds  stress  u'v'  ,  nondimensional ized  by  U^,  was  computed 
and  contour  plots  are  shown  in  Figures  28  through  31.  The  geometric 
pattern  is  qualitatively  similar  to  the  experimental  data  obtained  by 
Coles  and  Wadcock  (Reference  96)  for  a  NACA  4412  airfoil  which  are 
shown  in  Figure  32.  Some  experimental  measurements  were  made  during  the 
hot  wire  anemometry  study  (Reference  84)  of  the  velocity  profiles. 

This  data  also  illustrates  the  rapid  variation  of  negative  and  positive 
stresses  in  the  upper  and  lower  portions  of  the  near  wake,  respectively. 
Magnitudes  of  the  order  0.01  have  been  observed  (References  84,96)  in 
the  near  wake  and  were  also  obtained  numerically  where  the  contour  values 
vary  from  -0.01  to  +0.01  in  Figures  28  through  31. 

The  computed  surface  mean  pressure  distributions  for  all  four  angles 
of  attack,  presented  in  Figures  16  through  19,  favorably  compare  with 
the  experimental  data  (Reference  84)  which  have  an  estimated  error  of 
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5  to  8%.  An  inviscid  solution  with  imposed  Kutta  condition  was  computed 
using  the  approach  outlined  in  Appendix  G  and  is  also  shown.  The  computed 
Navier-Stokes  laminar  suction  peak  is  well  defined  and  increases  rapidly 
with  angle  of  attack.  The  lower  surface  pressure  peak  both  moves  down¬ 
stream  and  becomes  broader  as  the  angle  of  attack  increases.  This  behavior 
is  also  seen  in  the  experimental  data.  A  separation  bubble  is  observed 
in  the  experimental  data  where  the  curvature  rapidly  changes.  The  bubble 
occurs  further  downstream  when  compared  with  the  Navier-Stokes  computation. 
This  difference  is  probably  caused  by  freestream  turbulence  which  can 
delay  separation  and  the  subsequent  formation  of  the  bubble.  The  experi¬ 
mental  data  were  measured  in  a  wind  tunnel  with  a  freestream  turbulence 
intensity  of  approximately  0.5%  compared  with  the  unperturbed  freestream 
of  the  numerical  solution.  Some  disagreement  occurs  in  the  trailing  edge 
region  where  the  computed  result  has  a  more  pronounced  flat  trailing  edge 
stall  characteristic. 

The  computed  inviscid  surface  pressure  distribution  consistently 
predicts  both  a  larger  suction  peak  on  the  upper  surface  and  a  larger 
pressure  peak  on  the  lower  surface.  The  discrepancies  increase  with 
angle  of  attack.  The  inviscid  solutions  predict  complete  recovery  at  the 
trailing  edge  with  a  stagnation  point  and  =  1.0.  The  large  pressure 
gradients  in  the  inviscid  and  viscous  solutions  near  the  airfoil  nose 
are  resolved  which  indicates  an  adequate  grid  point  distribution  in 
the  region. 

An  examination  of  the  experimental  and  inviscid  surface  pressure 
distributions  for  zero  angle  of  attack  given  in  Figure  33  provides 
additional  information.  The  inviscid  pressure  distributions ,  obtained 
with  the  SOR  finite  difference  technique  in  Appendix  G,  for  the  upper 
and  lower  surfaces  are  indistinguishable.  This  physically  correct 
symmetric  result  came  directly  from  the  numerical  solution  and  was  not 
specified  as  a  boundary  condition.  Furthermore,  the  numerical  inviscid 
result  is  virtually  identical  with  an  inviscid  solution  (Reference  67) 
using  the  method  of  Theodorsen  (Reference  12).  The  experimental  data  for 
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both  surfaces  agree  well  with  the  inviscid  results  which  indicates  the 
small  effect  that  viscosity  has  on  the  pressure  distribution  for  a 
streamlined  symmetric  airfoil  at  zero  angle  of  attack.  This  agreement 
also  provides  evidence  that  the  experimental  airfoil  section  does 
approximate  the  NACA  0012  configuration.  Viscous  effects  on  the  pressure 
distribution  do  occur  (as  was  previously  mentioned)  as  the  angle  of  attack 
is  increased  even  for  small  a.  A  comparison  between  the  experimental 
data  for  the  two  surfaces  reveals  a  slight  suction  peak  on  the  lower 
surface.  The  suction  peak  suggests  that  the  experimental  zero  angle  of 
attack  may  actually  be  slightly  negative.  A  small  deviation  in  angle 
is  plausible  since  the  experimental  zero  was  determined  by  manually 
adjusting  the  airfoil  attitude  until  the  pressure  distributions  on  both 
surfaces  appeared  identical . 

The  lift  coefficients  obtained  from  the  present  turbulent  Navier- 
Stokes  solution,  inviscid  results,  and  two  sets  of  experimental  data  are 
compared  in  Figure  34.  The  experimental  results  shown  have  been  corrected 
for  angle  of  attack.  The  effective  experimental  zero  angle  of  attack  for 
a  symmetric  airfoil  section  is  defined  to  be  the  angle  where  the  lift 
coefficient  is  zero.  For  the  Seiler  Laboratory  results  (Reference  84), 
the  effective  zero  was  found  to  be  0.5°.  This  small  correction  is  in 
accord  with  the  analysis  of  the  nominal  zero  angle  of  attack  pressure 
distribution  results  already  presented.  Therefore,  the  experimental 
data  have  been  translated  0.5°  in  the  negative  a  direction  in  Figure  34. 
The  NACA  results  (Reference  56)  were  similarly  published  in  corrected 
form.  Other  common  two-dimensional  wind  tunnel  corrections  including 
solid  blocking,  wake  blocking,  and  streamline  curvature  effects  were 
investigated.  Empirical  low-speed  wind  tunnel  correction  factors 
(Reference  97)  using  the  Seiler  Laboratory  2'x3‘  tunnel  geometry  were 
applied  to  the  experimental  data  at  11.5°  to  account  for  blockage  and 
wall  interference.  The  corrected  values  for  the  lift  coefficient  and 
angle  of  attack  became  98%  of  the  uncorrected  lift  coefficient  and  +  0.1°, 
respectively.  Neither  data  set  given  in  Figure  34  have  been  so  corrected 
since  the  changes  are  small. 
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The  two  sets  of  experimental  data  contain  some  contrasts.  The 
NACA  data  (Reference  56)  maintain  an  approximately  linear  behavior  over  a 
greater  range  of  angles  of  attack  when  compared  with  the  more  recent  data 
(Reference  84).  The  NACA  lift  curve  then  flattens  more  rapidly  near  stall. 
These  minor  differences  may  be  explained  by  a  larger  freestream  turbulence 
level  in  the  NACA  results  which  tends  to  delay  the  turbulent  trailing 
edge  separation.  Also,  the  major  variation  of  the  lift  curve  at  small 
angles  of  attack  occurs  for  Reynolds  numbers  less  than  10^  (References 
56,89).  Freestream  turbulence  can  cause  more  rapid  transition  and  result 
in  an  apparent  freestream  Reynolds  number  greater  than  the  nominal 
turbulence-free  value. 

The  computed  turbulent  Navier-Stokes  lift  coefficients  for  the 
four  angle  of  attack  cases  are  in  excellent  agreement  with  the  experi¬ 
mental  data  as  shown  in  Figure  34.  The  force  components  were  computed 
using  trapezoidal  integration  to  sum  the  total  surface  stresses  obtained 
from  the  flow  field  solution.  The  numerical  results  exhibit  a  gradual 
decrease  in  slope  with  increased  angle  of  attack  similar  to  the  Seiler 
Laboratory  data.  This  behavior  is  consistent  with  the  data  because  the 
numerical  solution  has  a  quiescent  incoming  freestream  with  a  correspond¬ 
ing  smaller  apparent  Reynolds  number.  The  significant  effect  that  the 
viscous  separation  in  the  trailing  edge  region  has  on  the  lift  coefficient 
for  angles  of  attack  greater  than  8°  is  observed.  The  computed  values 
are  within  5%  of  the  experimental  results.  The  lift  coefficient  was  also 
computed  for  each  angle  of  attack  using  the  contour  momentum  integral 
method  described  in  Appendix  D.  The  calculated  values  differ  by  less 
than  1%  for  n  contour  paths  within  one-half  chord  of  the  airfoil  surface. 
This  result  demonstrates  the  near  flow  field  consistency  in  resolving  the 
lift  force.  The  time  dependent  terms  in  Equations  83  and  84  always  con¬ 
tributed  less  than  0.5%  of  the  total  computed  lift  which  indicates  again 
a  converged  steady  solution. 

Two  inviscid  flow  predictions  for  the  lift  coefficient  are  shown  in 
Figure  34  for  comparison.  The  numerical  inviscid  finite  difference  results 
were  computed  with  the  second  order  accurate  trapezoidal  integration 
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technique  used  in  the  viscous  calculations.  The  linear  airfoil  theory 
2 it  slope  is  included.  The  iuviscid  flow  results  overestimate  lift  by 
25%  at  an  angle  of  attack  of  5°.  This  discrepancy  increases  with  angle 
of  attack  because  the  viscous  effects  which  induce  stall  are  not  pre¬ 
sent  in  an  inviscid  calculation. 

The  computed  drag  coefficients  for  the  four  cases  are  compared  with 
available  experimental  data  (Reference  56)  in  Figure  35.  The  force  com¬ 
ponents  were  computed  using  trapezoidal  integration  to  sum  the  total 
surface  stresses  obtained  from  the  flow  field  solution.  The  drag  compo¬ 
nent  is  in  the  direction  of  the  incoming  freestream  flow.  The  rapid 
increase  in  the  drag  coefficient  which  accompanies  the  smaller  increase 
in  lift  at  higher  angles  of  attack  near  stall  is  observed.  Jacobs  and 
Sherman  (Reference  56)  acknowledged  that  the  experimental  drag  data  in 
this  Reynolds  number  region  have  error  greater  than  the  estimated  +0.001 
for  the  larger  Reynolds  number  results.  The  agreement  between  the  present 
numerical  solution  and  the  experimental  data  is  within  10  drag  counts 
in  the  region  of  the  maximum  lift  to  drag  ratio. 

The  presented  numerical  solutions  for  two-dimensional  incompressible 
turbulent  viscous  flow  over  airfoils  were  obtained  with  several  parametric 
studies.  The  effects  associated  with  varying  the  turbulence  transition 
values  and  the  SOR  iteration  parameters  including  time  step  size  have 
already  been  discussed.  The  influence  of  other  turbulence  model  para¬ 
meters,  far  field  boundary  conditions,  and  grid  fineness  on  the  numerical 
solutions  have  also  been  investigated. 

The  turbulence  model  contains  several  parameters  which  can  vary. 

The  effects  that  changes  in  these  parameters  have  on  the  flow  field  near 
the  airfoil  surface  and  on  the  values  of  lift  and  drag  were  examined. 

The  value  of  the  nominal  inner  eddy  viscosity  parameter  k-j  was  increased 
by  20%  at  a  5°  angle  of  attack.  The  lift  coefficient  subsequently 
increased  by  2%  while  the  drag  coefficient  increased  by  5%  or  six  drag 
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counts.  These  integrated  effects  are  caused  by  the  observed  small  increase 
of  velocity  in  the  inner  boundary-layer  on  the  upper  surface.  The  outer 
eddy  viscosity  parameter  kg  was  next  increased  by  50%.  The  result  was  an 
increase  in  turbulent  mixing  with  a  corresponding  thicker  boundary-layer. 
The  velocities  near  the  edge  of  the  upper  surface  boundary-layer  decreased 
which  caused  an  8%  loss  of  lift  with  no  apparent  effect  on  the  calculated 
drag. 


Values  for  the  constants  in  the  inner  eddy  viscosity  parameter  k1 
relaxation  factor  given  by  Equation  65  were  obtained  from  a  parametric 
study  at  an  angle  of  attack  of  11.5°.  This  angle  of  attack  was  chosen 
because  the  effects  of  the  relaxation  are  significant  only  for  angles  of 
attack  near  stall.  A  solution  was  initially  computed  without  the  relaxa¬ 
tion  factor  f.  The  laminar  separation  bubble  remained  closed  with 
trailing  edge  separation  beginning  at  mid-chord.  The  calculated  values 
for  the  lift  and  drag  coefficients  were  0.98  and  0.064,  respectively.  An 
examination  of  the  lift  curve  in  Figure  34  reveals  that  this  solution 
exhibits  the  character  of  leading  edge  stall  where  an  approximate  linear 
behavior  is  sustained  until  separation  abruptly  occurs.  Leading  edge  stall 
does  occur  for  the  NACA  0012  airfoil  for  Reynolds  numbers  greater  than 
about  500,000  (References  56,  89).  Thus  the  mechanism  for  leading  edge  stall  seems 
present  within  the  modified  turbulence  model.  A  numerical  solution  was 
next  obtained  using  the  relaxation  factor  with  a  relaxation  distance  sr  = 

0.25  and  a  delay  distance  s^  =  0.  in  nondimensional  distances.  The  laminar 
bubble  remain  closed,  but  the  trailing  edge  separation  region  was  signifi¬ 
cantly  larger  and  moved  to  within  the  20%  chord  location  (x  =  -  0.3).  The 
lift  coefficient  decreased  substantially  to  a  value  of  0.8  while  the  drag 
coefficient  increased  moderately  to  a  value  of  0.074.  The  final  solution 
which  has  been  previously  discussed  was  calculated  using  distances  $r  = 

0.5  and  s^  »  0.1.  The  computed  lift  and  drag  coefficients  became  0.84 
and  0.068,  respectively,  and  are  in  agreement  with  the  experimental  data. 

The  sensitivity  of  the  numerical  solution  to  these  parameters  has  thus 
been  obtained.  The  last  set  of  distance  parameters  was  used  in  calcu¬ 
lating  the  solutions  for  the  other  angles  of  attack.  The  onset  of 
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trailing  edge  separation  for  the  5°  and  7.5°  solutions,  seen  in  Figures 
4  and  5,  indicates  that  the  relaxation  factor  was  not  activated.  Thus, 
the  relaxation  model  for  parameter  affects  the  flow  field  only  near 
stall  where  the  pressure  gradient  approaches  zero  over  a  large  portion 
of  the  upper  surface. 

The  relaxation  distance  WLg  in  the  near  wake  turbulence  model  given 
by  Equation  74  was  also  investigated.  Numerical  solutions  were  obtained 
at  5°  angle  of  attack  for  relaxation  distances  equal  to  5  and  100  times 
the  lower  surface  boundary-layer  thickness  near  the  trailing  edge.  No 
changes  in  either  the  surface  pressure  distribution  or  the  integrated 
force  coefficients  were  observed.  In  two  recent  near  wake  calculations 
with  similar  relaxation  models,  Waskiewicz  (Reference  98)  used  a  value 
of  30  trailing  edge  boundary-layer  thicknesses  and  Hasen  (Reference  99) 
used  10  boundary-layer  thicknesses.  The  relaxation  distance  of  five 
boundary  layer  thicknesses  was  retained  since  the  agreement  of  the  Reynolds 
stress  field  discussed  previously  was  improved  for  this  choice. 

The  effect  on  the  solution  of  the  placement  and  type  of  far  field 
boundary  conditions  should  be  considered  in  any  computational  work.  In 
this  investigation,  two  types  of  boundary  conditions  at  different  loca¬ 
tions  were  examined.  The  four  angle  of  attack  solutions  that  have  been 
discussed  were  computed  using  the  freestream  far  field  boundary  conditions 
described  in  Sections  III. 2  and  IV. 3  at  a  circular  outer  boundary  of  radius 
10  chords.  The  effect  of  the  far  field  boundary  placement  on  an  inviscid 
solution  for  Joukowski  airfoils  was  previously  accomplished  (Reference  82). 
This  analysis  indicated  that  a  radius  of  10  chords  was  sufficient  as 
discussed  in  Section  IV. 3.  The  effect  of  the  boundary  placement  on 
the  present  viscous  solution  was  examined  by  using  an  alternate  outer 
boundary.  The  J  =  40  n  contour  with  a  semi-major  x  axis  of  4.77  and 
semi -mi nor  y  axis  of  4.64  was  chosen  as  the  new  outer  boundary  because 
it  approximates  a  circle  with  half  the  previous  radius.  The  numerical 
solution  was  then  obtained  for  a  =  5°  and  compared  with  the  solution 
using  the  10  chord  radius  outer  boundary.  The  average  values  for  the 
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lift  and  drag  coefficients  were  identical.  The  solution  with  the  closer 
far  field  boundary  had  a  small  oscillatory  behavior  with  variations  in 
the  lift  and  drag  coefficients  of  1%  and  2%,  respectively.  Further  ex¬ 
amination  of  the  near  surface  flow  field  revealed  that  the  laminar 
separation  bubble  had  a  small  scale  unsteadiness  which  locally  affected 
the  pressure  field. 

The  far  field  potential  flow  boundary  condition  model  developed  in 
Appendix  F  was  next  applied  with  the  outer  boundary  locations  of  10  chords 
and  about  five  chords  in  turn.  This  model  approximates  the  small  pertur¬ 
bations  from  the  freestream  conditions  for  the  velocity  and  pressure 
fields  at  a  large  but  finite  distance  caused  by  the  presence  of  a  body  in 
the  flow  field.  In  this  approach,  the  upstream  boundary  conditions  for 
both  velocity  and  pressure  became  the  calculated  values  from  the  inviscid 
potential  flow  model.  The  downstream  boundary  condition  for  pressure 
also  became  the  corrected  value  rather  than  the  freestream  value.  The  no 
change  downstream  boundary  condition  for  the  velocity  components  was  re¬ 
tained.  The  solution  was  calculated  with  the  revised  far  field  boundary 
conditions  at  the  10  chord  radius  circular  outer  boundary.  The  boundary 
conditions  were  initially  updated  every  200  time  steps  (0.2T)  by  using 
the  latest  calculated  value  of  the  circulation.  After  three  characteristic 
time  periods  T,  the  mean  values  of  the  lift  and  drag  coefficients  were 
virtually  identical  with  the  previous  freestream  boundary  condition  results. 
However,  while  the  lift  coefficient  had  variations  of  less  than  1%  from 
the  mean  value,  the  calculated  drag  coefficient  was  periodic  with  a 
period  of  0.4T  and  variations  of  +  10%.  Since  the  period  of  the 
oscillations  was  twice  the  period  of  the  boundary  condition  update 
and  the  circulation  was  essentially  a  constant  at  this  time,  the  outer 
boundary  conditions  (except  the  no  change  downstream  condition)  were  kept 
at  the  current  values  and  the  solution  was  advanced  1.5  characteristic 
time  periods.  At  this  time,  the  lift  coefficient  had  a  steady  value  of 
0.425  and  the  drag  coefficient  became  0.011  +  2%  (two  drag  counts). 

These  computed  results  are  within  1%  of  the  corresponding  coefficients 
calculated  from  the  solution  using  the  simple  freestream  conditions. 
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Thus  the  more  accurate  boundary  condition  which  includes  the  effect  of 
a  finite  distance  from  the  body  induces  very  small  changes  in  the  computed 
flow  field  near  the  airfoil. 

The  modified  far  field  potential  flow  boundary,  condi tons  were  also 
used  at  the  near-circular  five  chord  radius  outer  boundary  defined  previously 
to  obtain  a  solution  again  for  a  =  5°.  After  two  characteristic  times, 
the  outer  modified  freestream  boundary  conditions  were  held  constant  since 
the  lift  coefficient  had  small  variations  less  than  1%.  The  solution  ex¬ 
hibited  a  small  oscillatory  behavior  similar  to  the  freestream  boundary 
condition  result  even  after  four  characteristic  time  periods.  The  lift 
coefficient  had  a  mean  value  within  1%  of  the  freestream  boundary  condition 
value  with  variations  below  1%.  The  drag  coefficient  had  the  same  average 
value  as  the  previous  result,  but  the  oscillations  were  larger  with  a  de¬ 
viation  from  the  mean  of  approximately  +10%  (10  drag  counts).  The  un¬ 
steadiness  was  again  observed  to  emanate  from  the  laminar  separation  bubble. 
Small  pressure  oscillations  propagated  downstream  on  the  upper  surface 
of  the  airfoil  and  were  eventually  damped  out  near  the  trailing  edge.  In 
each  case  no  changes  were  made  in  any  of  the  turbulence  or  other  parameters 
during  the  computations.  Small  changes  in  the  turbulence  transition  near 
the  bubble  would  probably  prevent  the  unsteadiness. 

The  far  field  boundary  condition  study  has  indicated  that  the  use 
of  the  freestream  boundary  conditions  at  the  large  but  finite  distance 
from  the  airfoil  surface  is  sufficient  for  the  computation  of  the  near 
surface  flow  field  and  force  coefficients.  Variations  in  the  outer  bound¬ 
ary  placement  and  the  use  of  a  more  accurate  far  field  boundary  condition 
had  negligible  effects  on  the  present  numerical  solution. 

The  numerical  implementation  of  the  downstream  no  change  boundary 
condition  was  also  investigated.  Second  order  accurate  central  spatial 
and  upwind  differences  were  used  in  separate  solutions  for  an  angle  of 
attack  of  7.5°.  No  difference  (+0.01%)  was  observed  in  the  surface  pres¬ 
sure  distributions  or  computed  force  coefficients  between  the  two  bound¬ 
ary  conditions.  The  computed  velocities  at  locations  across  the  wake  on 
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the  outer  boundary  for  the  two  cases  differed  by  less  than  1%.  The  upwind 
difference  formulation  was  chosen  because  of  the  associated  convective 
properties. 

The  sensitivity  of  the  numerical  solution  to  the  fineness  of  the 
grid  was  examined.  The  7.5°  angle  of  attack  case  was  chosen  because 
both  the  laminar  separation  bubble  and  a  region  of  separated  flow  near 
the  trailing  edge  are  present.  A  coarse  grid  was  obtained  from  the  79x44 
grid  by  deleting  the  odd  numbered  n  contours  except  for  the  body  contour. 

The  coarse  79x23  grid  was  a  subset  of  the  full  coordinate  system  and  was 
used  to  study  the  effect  of  boundary-layer  resolution  on  the  numerical 
solution.  The  turbulence  parameters  and  SOR  acceleration  parameters  were 
unchanged  from  the  previous  solution.  The  coarse  grid  computation  was 
carried  out  over  four  characteristic  time  periods.  The  calculated  values 
for  the  lift  and  drag  coefficients  approached  a  "steady"  periodic  state 
after  two  characteristic  times  with  a  period  of  0.6T.  The  mean  value  of 
the  lift  coefficient,  averaged  over  the  last  two  characteristic  times, 
became  0.606  +  2%  compared  with  the  full  grid  solution  value  of  0.606. 

The  computed  mean  drag  coefficient  was  0.0260  +  15%  compared  with  the 
79x44  grid  solution  result  of  0.0240.  Values  for  the  wake  velocities  in 
the  freestream  direction  along  the  outer  boundary  for  each  solution 
compared  within  1%.  The  unsteadiness  was  found  to  originate  from  the 
oscillating  laminar  separating  shear  layer  which  forms  the  front  portion 
of  the  bubble.  The  edge  of  the  boundary-layer  varied  between  the  J=3 
and  J=4  n  contours,  as  defined  by  the  79x23  grid  system,  at  the  1%  chord 
location.  This  motion  also  caused  the  bubble  to  vary  in  length  and  produced 
pressure  oscillations  along  the  upper  surface.  This  unsteadiness  is 
probably  attributed  to  the  reduced  resolution  of  the  bubble  in  the  n 
direction.  Similar  sensitivity  has  already  been  discussed  concerning 
the  resolution  of  the  extremely  short  bubble  in  the  c  direction  at 
larger  angles  of  attack. 
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The  computed  flow  field  near  the  airfoil  and  the  calculation  of  the 
force  coefficients  using  the  coarse  grid  are  in  agreement  with  the  results 
for  the  complete  79x44  grid.  A  quadratic  Richardson's  extrapolation  on 
the  computed  drag  coefficients  indicates  an  error  of  +0.0007.  This 
agreement  indicates  that  the  solution  obtained  with  the  full  79x44 
coordinate  system  is  adequate  to  obtain  +10  drag  counts  and  therefore 
sufficiently  approximates  the  limiting  numerical  solution  obtained  from 
an  extremely  fine  grid.  The  coarse  grid  computation  further  indicates 
that  the  numerical  multi-grid  approach  (Reference  100)  may  be  implemented 
for  complex  flows  using  general  grid  transformations.  The  use  of  only 
one-half  the  grid  points,  for  instance,  during  most  of  the  time  marching 
procedure,  would  significantly  reduce  the  computer  time  required  for  an 
equivalent  numerical  solution. 
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SECTION  VI 

CONCLUSIONS  AND  RECOMMENDATIONS 

Numerical  solutions  have  been  obtained  for  two-dimensional  incom¬ 
pressible  turbulent  flow  over  airfoils  near  stall.  The  time  dependent 
Reynolds  averaged  incompressible  Navier-Stokes  equations  in  the  primitive 
variables  of  velocity  and  pressure  together  with  a  Poisson  pressure 
equation  are  numerically  solved.  An  algebraic  eddy  viscosity  approach 
is  modified  for  separated  adverse  pressure  gradient  flows  and  used  to 
model  turbulent  closure  of  the  laminar  separation  bubble  and  the  subse¬ 
quent  turbulent  boundary-layer .  A  deficiency  in  the  standard  model  is 
detected  and  corrected  by  using  a  "limiting"  technique.  A  body-fitted 
coordinate  system  is  numerically  transformed  to  a  rectangular  grid  in 
the  computational  plane.  The  set  of  transformed  partial  differential 
equations  is  solved  with  an  implicit  finite  difference  method.  Successive- 
over-relaxation  iteration  is  used  to  solve  the  system  of  linearized 
difference  equations  at  each  time  step. 

Numerical  solutions  are  presented  for  a  NACA  0012  airfoil  near  stall 
at  a  chord  Reynolds  number  of  170,000.  A  short  laminar  separation  bubble 
located  near  the  upper  surface  suction  pressure  peak  is  obtained.  Computed 
laminar  separation  bubble  characteristics  including  the  criteria  for 
separation,  bubble  type,  and  turbulent  transition  agree  with  empirical 
results.  Surface  mean  pressure  distributions  are  presented  and  found  to 
compare  favorably  with  experimental  data.  The  separation  bubble  is  ob¬ 
served  for  angles  of  attack  of  5°  and  7.5°.  For  larger  angles  of  attack, 
the  small  bubble  essentially  disappeared  within  the  numerical  resolution 
of  the  streamwise  grid  spacing.  The  steep  leading  edge  suction  pressure 
peak  is  well  defined  for  each  angle  of  attack.  Velocity  profiles  at  four 
stations  along  the  upper  airfoil  surface  are  compared  with  experimental 
results.  The  experimental  data  were  obtained  at  similar  grid  point 
locations  so  that  interpolation  was  not  required.  The  calculated  lift 
and  drag  coefficients  are  in  excellent  agreement  with  the  experimental 
data.  The  lift  coefficients  are  within  5%  of  the  experimental  values 
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near  stall,  and  the  computed  drag  coefficients  are  within  ten  drag  counts 
in  the  region  of  maximum  lift  to  drag  ratio.  The  effects  of  viscous 
separation  on  the  lift  coefficient  curve  which  produces  a  maximum  value 
are  seen  and  cjntrasted  with  a  numerically  obtained  inviscid  potential 
solution  with  Kuttn  condition.  The  observed  phenomena  of  trailing  edge 
stall  is  predicted  where  the  rear  separation  point  moves  forward  with 
increasing  angle  of  attack. 

The  sensitivity  of  the  numerical  solution  has  been  examined  in 
several  areas.  A  far  field  potential  flow  boundary  condition  which 
modifies  the  freest ream  conditions  for  use  at  a  large  but  finite 
distance  from  the  airfoil  was  considered  for  two  outer  boundary  place¬ 
ments.  The  study  showed  that  the  use  of  the  infinite  freestream  condi¬ 
tions  at  an  outer  circular  boundary  of  radius  ten  chords  produces 
negligible  influence  on  the  near  flow  field  and  calculated  values  for 
the  force  coefficients.  A  coarse  79x23  grid,  compared  with  the  79x44 
grid,  was  used  to  evaluate  truncation  error.  An  analysis  of  convergence 
for  successive-over-relaxation  iteration  predicts  an  upper  limit  on  the 
time  increment  for  the  implicit  finite  difference  method.  Numerical 
experiments  confirmed  this  upper  bound.  Several  parameters  within  the 
turbulence  model  were  systematically  varied.  The  only  significant 
sensitivity  occurred  near  the  downstream  side  of  the  laminar  separation 
bubble.  Small  changes  in  turbulent  transition  length  and  location  were 
found  to  significantly  change  the  flow  field.  This  sensitivity  may  be 
related  to  incipient  turbulent  separation  which  results  in  failure  to 
close  the  separation  bubble.  This  leading  edge  stall  phenomena  is 
observed  at  higher  Reynolds  numbers  for  the  NACA  0012  airfoil. 

The  agreement  between  the  numerical  solutions. and  the  experimental 
data  and  empirical  results  indicates  that  the  eddy  viscosity  approach 
modified  for  separated  adverse  pressure  gradient  flows  adequately 
models  the  turbulence  on  the  upper  airfoil  surface.  The  near  stall 
airfoil  aerodynamic  characteristics  are  consequently  accurately  pre¬ 
dicted  by  the  numerical  method. 
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The  results  of  this  investigation  have  suggested  areas  for  further 
research.  A  more  exhaustive  study  of  the  laminar  separation  bubble 
should  be  accomplished  to  better  understand  the  phenomena  of  leading  edge 
stall.  The  computation  of  a  geometry  with  detailed  experimental  results 
for  this  purpose  is  desirable.  Further  research  in  adapting  the  multi¬ 
grid  technique  would  significantly  increase  the  computational  efficiency 
and  appears  feasible  as  a  result  of  the  coarse  grid  calculation. 
Turbulence  modelling  is  an  area  of  substantial  interest  and  should  be 
pursued  in  conjunction  with  experimental  investigations.  The  improved 
accuracy  of  current  computational  aerodynamics  also  indicates  that  the 
errors  in  experimental  measurements  should  be  analyzed  in  a  more  rigorous 
manner  so  that  comparisons  can  be  better  evaluated.  The  method  should 
also  be  employed  in  time  dependent  problems  to  further  exploit  the  time 
marching  technique. 
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Figure  1.  Physical  and  Transformed  Planes 
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Figure  3.  Portion  of  79x44  Body-Fitted  Grid  (First  30 
Constant  ri  Lines) 
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Figure  7.  Mean  Flow  Streamline  Contours  ( 
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Figure  8.  Mean  Flow  Velocity  Field  (i 


Figure  9.  Mean  Flow  Velocity  Field  ( 


Figure  10.  Mean  Flow  Velocity  Field  ( 
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Figure  11.  Mean  Flow  Velocity  Field  (a 


Figure  12.  Mean  Flow  Pressure  Contours  ( 


Figure  13.  Mean  Flow  Pressure  Contours  (a 
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Figure  14.  Mean  Flow  Pressure  Contours  (a  =  9.5°) 


Figure  15.  Mean  Flow  Pressure  Contours  (a  =  11.5°) 
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Figure  19.  Surface  Mean  Pressure  Coefficients  ( 
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Figure  20.  Mean  Flow  Velocity  Versus  Surface  Normal  Distance 

Nondimens ional ized  by  Computed  Boundary-Layer  Thick¬ 
ness  at  Chord  Locations  on  Upper  Surface  (a  =  5°) 
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Figure  21.  Mean  Flow  Velocity  Versus  Surface  Normal  Distance 

Nondimensional ized  by  Computed  Boundary-Layer  Thick¬ 
ness  at  Chord  Locations  on  Upper  Surface  (a  =  7.5°) 
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Figure  22.  Wean  Flow  Velocity  Versus  Surface  Normal  Distance 

Nondimensionalized  by  Computed  Boundary-Layer  Thick, 
ness  at  Chord  Locations  on  Upper  Surface  (a  =  9.5°) 


Figure  23.  Mean  Flow  Velocity  Versus  Surface  Normal  Distance 

Nondimensionalized  by  Computed  Boundary-Layer  Thick¬ 
ness  at  Chord  Locations  on  Upper  Surface  (a  =  11.5°) 
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Figure  26. 


Wake  Mean  Flow  Velocity  Profiles  at  Chord  Locations 
with  Origin  on  Chordline  (a  =  9.5°) 


Figure  27.  Wake  Mean  Flow  Velocity  Profiles  at  Chord  Locations 
with  Origin  on  Chordline  (a  =  11.5°) 
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Figure  30.  Reynolds  Stress  o'”'  rnntours  (a  =  9.5.°) 
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Figure  31.  Reynolds  Stress  u’v’  Contours  (a  *  11.5°) 
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Figure  II.  Contour  plot*  of  double  correlation*  <o'u'>/(q 
<u*v'>/(qr€f)*,  »nd  <vlv'>/(qre^)*.  Contour  interval  0. 010  or  0. 005. 


Figure  32. 


Experimental 
4412  Airfoil 


Reynolds  Stress  u'v'  Contours  for  NACA 
(a  =  14°  and  Re  =  1 .5  x  10^),  Ref.  96 


Figure  33.  Surface  Mean  Pressure  Coefficients  Comparing  Numerical 
Inviscid  Methods  with  Experimental  Data  at  a  =  0° 


108 


AFVIAL-TR-81 -3053 


Figure  34.  Lift  Coefficient  Curve 
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Figure  35.  Drag  Polar  Curve 
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CROSS-SECTION  AREA  A 


Figure  36.  Contour  Integration  Geometry  for 
Two-Dimensional  Flow 
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Figure  37.  Maximum  Error  of  SOR  Solution  for  U,  P,  P^  Versus 
Number  of  Iterations  with  at  =  0.0005  (a  =  8°) 
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Figure  38.  Maximum  Error  of  SOR  Solution  for  U,  P,  P^  Versus 
Number  of  Iterations  with  At  =  0.001  (a  =  8°) 
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APPENDIX  A 

GRID  TRANSFORMATION  RELATIONSHIPS 


Several  of  the  definitions  and  coordinate  transformation  relations 
in  the  transformed  plane  are  summarized  below.  The  notation  used  by 
Thames  (Reference  44)  and  Hodge  (Reference  53)  is  retained.  The 

physical  plane  is  the  (x,  y)  coordinate  system  and  the  computational 

plane  is  the  (c,  n)  coordinate  system.  The  following  function 
definitions  are  used: 

f(x,  y,  t)  -  a  twice  differentiable  scalar  function. 

F(x,  y,  t)  =  i  F-j(x,  y,  t)  +  j  F2  (x,  y,  t)  a  twice  differentiable 

vector  function  where  i  and  j  are  Cartesian  unit  vectors. 

Transformation  Definitions 


2  2 

0  =  x  +  y 
n  Jn 

(A-l) 

8  =  Vn  +  Vn 

(A-2) 

2  2 
y  =  xc  +  yc 

(A-3) 

D1  =  Kc  '  28  X5n  +  yXnn} 

(A-4) 

°2  =  Kc  '  28  y?n  +  Yynn} 

(A-5) 

J  =  x?yn  -  Xny? 

(A-6) 

°  = 

(A-7) 

T  =  <xr,°2  -  yn°l)/J 

(A-8) 

Derivatives  of  Scalar  Functions  in  Transformed  Plane 


fx  -  af/a*  -  (y/5  -  Vn>/J 


(A-9) 
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fy  =  af/ay  =  (x^  -  x^/J  (A-10) 

f.  =  3f/3t  =  (3f/at)  _  for  fixed  coordinates  (A-ll) 

t  S,n 

Vector  Operators  in  the  Transformed  Plane 
Gradient: 

-f  ■  -  *ev 5  +  <v„  -  \V  h/J 

Divergence: 

I  •  £  *  iy„(F,)  5  -  ye(F,l  „  *  xe(f2)  n  -  Xn(F2)  £)/J  (fl-13) 

Curl : 

:>i=  My„(F2>  E  -  Vf2)  n  -  x5(F, )  n  4  xn(F, )  51/J  (fl-14) 

Laplacian: 

v2f  •  («f?5  -  26 f  +  YfnT))/J2  +  (a/J2)fn  +  (x/J2)fe  (A-15) 

Unit  Normal  and  Tangent  Vectors  in  Physical  Plane 
Normal  to  n  line: 

=  vn/|vn|  =  (-y^i  +  x^j/ZT  (A-16) 

Normal  to  c  line: 

=  vc/|vc|  =  (y^i  -  x^j)//-^  (A-17) 

Tangent  to  n  line  (increasing  £  direction): 

^(n)  _  ^(n)  x  £  _  (x^  +  y^ j )/ i/~~y  (A-18) 

Tangent  to  c  line  (increasing  n  direction): 

_t^  =  -n^  x  k  =  (x^i  +  yrij)//”ci  (A-19) 
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Integrals  in  the  Transformed  Plane 
Area  Integral : 

/Rf(x.  y)  dxdy  =  /R*  f(x(?,  n),  yU,  n))  |J|  d?dn 

where  R  Is  the  region  R  mapped  into  the  (?,  n)  plane. 
Vector  Line  Integral: 


/^max 

F(x(5, 

C  . 
min 


op)>  y(C»  n  ))/"y  d£ 


(A-20) 


(A-21 ) 


where  contour  s^  is  any  constant  n  line  in  the  physical  plane,  s  is 
arclength  along  s^ ,  and  np  is  the  value  of  n  for  the  chosen  contour 
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APPENDIX  B 

FINITE  DIFFERENCE  APPROXIMATIONS 

This  appendix  presents  the  finite  difference  approximations  that  are 
used  in  this  investigation.  The  differences  are  formulated  using  a  func¬ 
tion  fU,  n)  in  the  computational  or  transformed  plane.  The  A?  and  An 
spacings  are  assumed  constant  with  value  unity.  The  truncation  error  term 
is  given  in  derivative  form  assuming  a£  and  An  are  equal.  The  approxima¬ 
tions  are  second  order  accurate  unless  otherwise  specified.  Time  dif- 

ferences  are  expressed  by  At.  The  superscript  n  denotes  the  nL  time 
interval  and  is  understood  when  omitted.  The  differences  are  given  for  a 
(£,  n)  point  location  denoted  by  subscripts  (i,  j)  and  are  understood  when 
omitted.  Space  derivatives  with  respect  to  only  n  are  presented  because 
the  corresponding  5  derivatives  are  identical  with  subscripts  (i,  j) 
reversed. 


First  time  derivative,  backward  difference,  f i-st  order: 
ft  =  (fn  .  fn_1)/At  -  At  fu/2  +  ... 


First  derivative,  backward  difference,  first  order: 

f  =  (f.  -  f.  .)  -  f  /2  +  ... 
n  j  j-1  nn 

First  derivative,  central  difference: 

f  =  '  f-  7)/2  -  f  /6  +  ... 

n  j+1  J-'  nnn 


First  derivative,  forward  difference: 

f  =  (-f +  4f.Al  -  3f.)/2  -  f  /3  +  ... 
n  j+2  j+1  j  nnn 


First  derivative,  backward  difference: 

f  =  (f <  ?  -  4f  .  +  3f.)/2  +  f  /3  +  ... 
n  i-<-  j-i  j  nnn 


(B-l ) 

(B-2) 

(B-3) 

(B-4) 

(B-5) 
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Second  derivative,  central  difference: 

f  =  (f...  -  2f.  +  f.  ,)  -  f  /1 2  +  ... 
nn  j+1  j  j-1  nnnrj 

Second  derivative,  forward  difference: 

f  ■  (-f,-xo  +  4f,x,  -5f .  .  +  2f.)  +  Ilf  /1 2  +  ... 
nn  j+3  j+2  j+1  j'  nnnn 

Second  derivative,  backward  difference: 

nn  j-3  j-2  j-1  y  nnnn 

Cross  derivative,  central  difference: 

f5n  "  (f1+l,  j+1  "  f1+1,  j-1  *  f1-l,  j-1  ‘  fi-l,  j+l)/4 

-  (f«En  +  W/24  *  - 

Cross  derivative,  central  in  £  and  forward  in  n  differences: 


(B-6) 


(B-7 ) 


(B-8) 


(B-9) 


j+i  ■  fi-i,  j+i’ 


■  <-3(ft+i  -  fi-,)  t4<fi+i, 

-  (f1+l.  j+2  -  ft-l,  j+2),/4  +  (fEEE„  *  W/6+' 


(B-10) 


Cross  derivative,  central  in  £  and  backward  in  n  differences: 
fEn  ‘  13<fi+l  -  f1-1>  -  «<',♦!.  j-1  -  f,-l.  j-,) 

t(fi+l,j-2-fi-l.j-2>l24-(fE«n*fEnnn>/6* 


(B-ll) 


117 


AFWAL-TR-81  -3053 


APPENDIX  C 

GRID  MODIFICATION  FOR  WAKE  RESOLUTION 

This  appendix  describes  the  technique  that  is  used  to  insert  addi¬ 
tional  grid  lines  in  the  region  of  the  airfoil  wake.  The  body-fitted  coor¬ 
dinate  grid  system  becomes  coarse  at  distances  far  from  the  body.  The 
coarse  grid  may  reduce  the  resolution  of  the  velocity  defect  in  the  viscous 
far  wake  region.  To  increase  wake  resolution,  eight  additional  %  lines 
are  placed  downstream  of  the  airfoil  with  body  points  on  the  rounded 
trailing  edge.  The  original  numerically  generated  grid  has  71  5  and  44  n 
lines.  Linear  interpolation  between  adjoining  ?  lines  of  the  form  x  (I,J) 

=  { x ( I -1 ,  J)  +  x ( I ,  J)}/2  is  used  to  locite  a  C  line  between  lines  1=1 
and  I  =  2,  I  *  2  and  I  =  3,  I  =  69  and  I  =  70,  and  I  =  70  and  I  =  71 .  The 
procedure  is  repeated  using  the  new  lines  as  well  except  only  one  addi¬ 
tional  line  is  located  below  the  I  =  1  cut.  In  this  way  the  fine  grid 
spacing  has  an  asymmetry  for  better  wake  resolution  at  angle  of  attack. 

The  final  79  by  44  grid  is  shown  in  Figure  3.  Table  C-l  gives  the 
numbered  I  index  (5  line)  designation  for  the  grid  systems  in  the 
wake  region. 
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Table  C-l 

fc 

i 

I  INDEX  U  LINE) 

DESIGNATION  FOR 

GRID  SYSTEMS 

IN  MAKE 

■  ) 

r 

71  x  44  System 

75  x  44  System 

79 

x  44  System 

69 

71 

72 

tv 

- 

72 

73 

k; 

- 

- 

74 

70 

73 

75 

£ 

- 

- 

76 

- 

74 

77 

M  f 

- 

- 

78 

71,  1 

75,  1 

79,  1 

- 

- 

2 

r 

- 

2 

3 

2 

3 

4 

- 

4 

5 

3 

5 

6 
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APPENDIX  D 

CONTOUR  INTEGRAL  METHOD  FOR  DETERMINING  BODY  FORCE  COEFFICIENTS 


The  determination  of  the  resultant  force  that  a  flow  field  exerts  on 
a  body  usually  involves  the  calculation  of  surface  pressure  and  viscous 
stresses  and  the  summation  of  these  forces  over  the  body  surface.  An 
alternate  approach  is  to  apply  a  control  volume  analysis  to  a  region 
enclosing  the  body.  In  the  foregoing  discussion,  a  control  volume  analysis 
for  a  body  in  a  two-dimensional  flow  is  presented. 


Consider  a  general  body  immersed  in  a  flow  field  with  surface  Og. 

Define  a  fixed  control  volume  V  with  outer  surface  a  and  inner  surface  aD 

n  B 

Conservation  of  linear  momentum  for  the  fluid  in  V  is  expressed  by  the 
Cauchy  Integral  Equation  of  Motion 


D_ 

Dt 


pb  dV  + 


/ 


n  ‘  T  do 


( D— 1 ) 


a 

where  p  is  the  density,  v^  velocity  field,  b  body  force  per  unit  mass,  n 

unit  outward  surface  normal  vector,  T  traction  stress,  and  o  is  the  total 
surface  Og  +  o^.  Apply  the  Reynolds  Transport  Theorem  to  the  material 

derivative  term  and  the  Divergence  Theorem  to  the  surface  integral  and  obtain 

from  Equation  D-l  for  the  i^  Cartesian  vector  component 

y*|T(pv.)dV+  f  d.{  pv.  Vj  -  T^)  dV  =  0  (D-2) 

V  V 


Apply  the  Divergence  Theorem  to  the  second  volume  integral  and  express  the 

resulting  surface  integral  as  separate  integrals  over  o0  and  a  to  obtain 

d  n 
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f  h  (pvi)  dv  + 

V 


/  "j<pv1  *J  -  V  d°  ♦  f  n.jHvj 

“8  °n 


Vd° 


=0 


(D-3) 


Now,  on  the  solid  nonporous  body  surface  oD,  n.  v.  =  0.  The  integral 

d  j  j 

of  the  resultant  traction  force  over  ag  is  the  net  force  exerted  by  the 

body  on  the  control  volume  fluid.  Therefore,  the  force  £  exerted  by  the 
fluid  on  the  body  is  given  by 

F,  •  -  /  n,  T.j  do  (0-4) 

Substitute  Equation  D-4  into  D-3  and  use  the  above  surface  condition  to 
obtain  a  general  expression  for  the  force  exerted  by  the  fluid  on  the  body. 


Fi  = 


at 


J"  PV|  dV  +  f 

I#  •/ 


nj(Tji  -  »Vj>  d0 


(0-5) 


Now  consider  the  case  of  two-dimensional  flow  in  the  x-y  plane.  Then 
Equation  D-5  becomes  the  following  expression  for  the  force  per  length  of 
span  £ 

fi  ■  -  ft  /  cVi dA  ♦  /  nj|Tjf  -  pVj)  ds  (0-6) 

n 


where  A  is  the  cross-sectional  irea  of  the  control  volume  and  s  is  the  outer 

n 

perimeter  of  A  as  shown  in  Figure  36. 


Substitute  for  the  traction  stress  T^  =  "P6ji  +  Tji  *n  Equation  D-6 


and  obtain 


f,  •  /  *  'j,  -  ‘>Vj)ds  -  ft  /  p,i  dA  (D-7) 
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where  p  is  the  pressure,  -r^  is  the  laminar  viscous  stress  and  <5^ 

is  the  Kronecker  Delta  function.  Next  express  the  flow  variables  in 
terms  of  turbulent  mean  and  fluctuating  variables,  assume  incompressible 
flow,  and  Reynolds  average  Equation  D-7  to  obtain 


/  njC-PJj,  *  xJf 
n 

J  pVj  dA 


-  p(vivj  +  vivj)}  ds 


(D-8) 


The  term  -pv'.v is  the  turbulent  Reynolds  stress  t.  ...  Let  the  total  stress 
■  J  ^  J  * 


be  defined  as  x^.  =  +  x^..  Then  Equation  D-8  becomes 

T,  -  /  n.  {-pjj,  *  tTjj  -  pVj}  ds  -|t  / 


dA  (D-9) 


Introduce  the  following  nondimensional  variables  in  Equation  D-9 


P  -  P„ 

- 2 — 


tu 


Vi  -  X1 

p  ./■  v.  =  jt- —  x.  =  t-L  and  t  =  -r—  ,  and  note  that 

pi)  1  U  T  L  L 


x,..  =77  . .  where  U  is  the  freestream  velocity  and  L  is  a  reference 

Tji  Tji  ®> 

f  •  i 

length  in  the  x  direction.  Also,  define  the  orce  coefficient  C.  *  X* 

•  2 


1  Pol¬ 


and  obtain  from  Equation  D-9 
C 


f-  2  *  Sin  ttji  •  Vj>  ds  - 2  Tt  /* idS  (,>-,0, 


where  all  variables  are  understood  to  be  mean  dimensionless  variables.  Now 
the  turbulent  eddy  viscosity  concept  Is  used  where 

xTji  *  +  eM^  *  ;  1vj  +  3jv1) 


and  define  Re 


PU  L 


t  u  + 


which  transforms  Equation  D-10  to  the  following 
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/ 


*  Lt  <yj 


+  Y<’ 


v.v.}  ds 

*  J 


2  It  /'idA 

A 

(D-ll) 


Expand  Equation  D-11  into  x  and  y  component  equations 

V  2  /  ('"'P  *  [2nlH  *  "2  <57  +  -  (“2"1  *  "Vn2)}  dS 

n 

-  2  |-t  J  udA  (D-12a) 

A 

Cfy*  2  f  {-n2p  +  ^  [n]  (fy  +  f^)  +  2n2  ”  (uvnl  +  v2n2)}  ds 

"  2  It  /  vdA  (0-1 2b) 

A 

where  n1  and  n2  are  the  x  and  y  components  of  the  unit  outward  normal 

a 

vector  n  for  the  outer  boundary  s^  of  area  A.  Equations  (D-12a)  and  (D-12b) 

are  expressions  for  the  x  and  y  force  coefficients  exerted  by  a  two- 
dimensional  incompressible  turbulent  flow  on  an  arbitrary  two-dimensional 
body  where  path  s^  encloses  the  body  and  an  eddy  viscosity  method  is  used 

to  model  turbulence. 


Lift  and  drag  coefficients  are  then  obtained  from  the  force  coefficients 
as  follows  where  a  is  the  geometric  angle  of  attack. 


C.  =  C*  cos  0  -  Cf  sin  a 
y  x 

(D-13) 

Cp  =  C-  Sin  a  +  C,  cos  a 

Ty  Tx 

(0-14) 
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APPENDIX  E 

TIME  INCREMENT  INFLUENCE  ON  SOR  CONVERGENCE 

Implicit  SOR  methods  have  been  shown  by  linear  stability  analyses 
to  be  unconditionally  stable.  However,  convergence  has  required  a  small 
time  step  size  in  this  investigation  and  elsewhere.  The  following  analy¬ 
sis  indicates  the  relationship  between  the  time  step  size  and  convergence. 

Consider  the  system  of  difference  equations  to  be  written  in  the 
following  matrix  form: 


AU  =  b 


(E-l) 


where  A  is  the  coefficient  matrix,  U  is  the  vector  of  unknown  primitive 
variables  at  time  step  n  with  u,  v,  p  grouped  together  for  each  (i,j) 
point  location,  and  b  is  a  vector  of  "constants".  For  example,  the  x  com 
ponent  of  the  momentum  equation  (Equation  92)  at  the  (i,j)  point  is 


q  u(*l  =  U?;1  ♦  att-YETA  <p)+,  j  -  P,.,  j)  *  YXI  (p,  j+,  -  p,  J.,) 

+  | UC|  (4uicl  j  -  UIC2  j)  +  I  VC  I  (4ui  JC1  -  u.  JC2)  +  |UV|  (4uJvl  . 

UIV2  +  lVVl  (4ui  JV1  *  ui  JV2^  +  RlT  {ALFA  (ui+l  j  +  ui-l  j* 


+  GAMA  (u,  j+1  +  u.  ._,)  +  BETA  (u.+]  j+]  -  u,+1  ♦  u..,  - 

«  (e-2) 


where  all  coefficients  are  evaluated  at  point  (i,j)  and  primitive  variables 
are  at  time  step  n  unless  given  otherwise,  and  Q  is  given  by 

Q  =  1  +  3  &t(|UC|  +  |VC|  +  |UV|  +  JV»|)  +  ( | ALFA |  +  J  GAMA | ) 
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Rearrange  Equation  E-2  and  obtain 

-  |UC|  Uic2  .  -  I UV |  U]v2  J  *  Ilf  uf_,  J.,  ♦  4  |UC|  U)cl  J  *  4  |U»|  UjV]  J 


“1-1  j  -  VETA  p,,,  j  -  jjff  Vl  J4l  -  'VCI  “1  OCZ  -|VV'  "1 


♦fSf  “1  1-1  "  l«l  “1  JC1  “1 


YXI  p.  .  .  -  J  u(.s) 
*1  J-I  At  13 


.  GAMA  „  .  V¥T  n  BETA  „  ,  ALFA  „  vrT.  , 

RET  Ui  j+1  YXI  Pi  j+1  '  RET  ui+l  j-1  RET  Ui+1  j  "  YETA  pi+l  j 

+  MAU  .  .  u"'1*  Js-1)  D  (E.3- 

RET  ui+l  j+1  rrji  +  Uij  D 


A  similar  y  component  of  momentum  equation  is  obtained  if  u  is  replaced  by 
v;  -  YETA  by  XETA;  and  YXI  by  -  XXI. 

The  finite  difference  pressure  Equation  94  becomes 

-  1 UV |  PIV2  j  ♦  BETA  p,_,  *4|UV|  p,#,  }  *  ALFA  p,_,  j 

-  B™  “1-1  jn  '  !vvl  “1  jV2  *  Pi  j-i+4lvvl  “I  jvi  -  % 

1  GAMA  p,  jt,  -  BETA  pH,  ,  4  ALFA  p,,,  }  *  BETA  p,+,  ,  -  -RHS 


*  O  p 

where  B=2  (ALFA+GAMA)  +  3  (|UV|  +  jw|)  and  RHS  =  -  ^-r-  +  C( UX ) ^+2 ( VX ) ( UY )+( VY )  ] 


Convergence  theorems  (Reference  81)  for  systems  of  linear  equations 
with  constant  coefficients  show  that  a  necessary  and  sufficient  condition 
for  convergence  requires  that  the  coefficient  matrix  be  diagonally  dominant. 
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If  this  criteria  is  applied  to  the  above  quasi-1 inear  system,  the 
following  convergence  criteria,  inequality  E-6  from  the  momentum 
equation  and  inequality  E-8  from  the  pressure  equation,  are  obtained. 

5  (|UC|  +  |VC|  +  |UV|  +  |VV|)  +  R|T  { |ALFAj+|GAMA| ) 

+  |ET  |  BETA  J  +  2(|YXI|  +  |  YETA | )  <  £t 

or  <  3  (|UC|  +  |VC|  +  |  UV|  +  |W|)  +  R|T  ( |  AL  FA  |  + 1  GAMA  | )  +  ^  (E-5) 

which  becomes 

•5  -1 

At<{2  (|UC|+  |VC|  +  |UV|  +  |VV|)  +  RfT|BETA|  +2( | YXI |+| YETA | )  (E-6) 

and 

2(  |  ALFA |  +  |  GAMA | )  +4  | BETA |  +5  (jUV|  +  |W|)  < 

2( | ALFA |  +  | GAMA | )  +  3  (|UV|  +  1 VV | )  (E-7) 

or 

4  | BETA |  +  2  (|UV]+|VV| )  <  0  (E-8) 


The  convergence  criteria  (Equation  E-6)  provides  a  limit  on  At  as  a 
function  of  both  the  grid  geometry  and  flow  field  solution.  The  right-hand 
side  of  the  inequality  has  its  smallest  magnitude  near  the  trailing  edge 
and  increases  rapidly  away  from  the  airfoil.  An  order  of  magnitude  analysis 
of  the  grid  coefficients  and  an  examination  of  the  computational  solutions 
Indicate  that  diagonal  dominance  occurs  if  At  <  0.0003  for  points  near  the 
trailing  edge.  At  <  0.001  elsewhere  near  the  body,  and  At  <  1  In  the  far 
field.  The  pressure  equation  convergence  criteria,  (Equation  E-8)  has  no 
At  dependence  and  Is  approximately  satisfied  only  at  points  far  from  the 
body  where  the  left-hand  side  approaches  zero.  These  criteria  can  only 
Indicate  a  possible  time  step  restriction  since  the  system  of  equations 
has  variable  coefficients  and  lower  order  nonlinearities. 
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A  numerical  experiment  was  conducted  to  determine  the  time  step  effect 
on  convergence.  A  computation  for  the  flow  over  the  NACA  0012  airfoil  at 
an  angle  of  attack  of  eight  degrees  was  used.  The  solution  was  advanced 
one  time  step  with  the  time  step  size  and  number  of  iterations  per  time 
step  varied.  The  solution  converged  for  time  steps  of  0.0005  and  0.001 
while  diverging  for  a  time  step  of  0.005.  The  relative  maximum  error 
magnitudes  for  the  u  component  of  velocity,  surface  pressure  p^,  and  field 
pressure  p  as  functions  of  number  of  iterations  are  shown  for  each  time 
step  in  Figures  37,  38,  and  39.  The  relative  error  is  defined  as 
(f  -  f  i  )/f  .  The  v  component  of  velocity  errors  behave  similarly 

to  the  u  component.  Next,  the  pressure  was  held  fixed  and  the  convergence 
of  the  velocity  was  monitored  {Figures  37,  38,  and  39).  The  rapid  rate  of 
convergence  is  observed.  The  u  and  v  velocity  components  were  also  held 
fixed  in  turn  to  test  convergence.  In  both  cases,  convergence  was  similar 
to  the  general  solution  case.  These  numerical  experiments  indicate  that 
the  convergence  of  pressure  is  slow  and  that  convergence  is  dependent 
on  time  step  sizes  similar  to  those  predicted  by  the  simple  analysis 
above. 


AFWAL-TR-81  -3053 


APPENDIX  F 

FAR  FIELD  BOUNDARY  CONDITIONS 


The  far  field  boundary  for  a  numerical  solution  must  usually 
be  at  a  finite  distance  from  the  body  of  interest.  This  constraint 
may  pose  difficulties  if  the  only  known  far  field  boundary  conditions 
are  those  at  infinity.  This  analysis  uses  complex  variable  methods 
for  incompressible  potential  flow  and  develops  approximate  far  field 
boundary  conditions  for  use  at  a  finite  distance. 


Consider  a  circular  cylinder  radius  a  and  center  at  the  origin,  with 
positive  clockwise  circulation  r,  in  a  uniform  stream  (velocity  U  ), 

oo 

at  angle  of  attack  a.  This  coordinate  system  is  identical  with  the 
physical  plane  of  the  airfoil.  The  complex  velocity  is 


W  (Z)  =  U  e"iot 

oo 


(2Tra2U  )  ia 

oo  0 

2ir  ^T 


+ 


( F-l ) 


where  the  first  term  is  the  uniform  stream  contribution,  the  second 
term  is  a  doublet  of  strength  2ira2Uco,  and  the  third  term  is  a  vortex  of 
strength  r. 


Next,  use  the  Kutta-Joukowski  Theorem  for  lift  and  the  definition 
of  the  lift  coefficient  to  obtain 

r  =  (1/2)  C^c  (F-2) 

where  c  is  the  airfoil  chord. 


Nondimenslonal ize  Equation  F-l  velocities  with  the  freestream  velocity 

and  lengths  with  the  airfoil  chord  c  and  substitute  Equation  F-2  into  F-l 
and  obtain 
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W  (Z)  =  e"1a 


_k  I 

4tt  Z 


( F-3) 


where  W  is  the  nondimensional  complex  velocity  and  Z  is  the  nondimensional 
complex  variable  Z  =  x  +  iy. 


Airfoil  transformations  such  as  the  Joukowski  and  Treffetz  Trans¬ 
formations  transform  a  cylinder  into  an  airfoil  with  a  chord  of  four 
radii.  Thus,  for  a  given  chord  c,  choose  the  cylinder  radius  to  be  1/4 

c  which  gives  a  scaling  factor  ~  =  1/4  in  the  doublet.  Then  the  complex 

velocity  for  the  flow  around  the  scaled  cylinder  becomes 


W  (Z)  =  e'lct 


Jot 


W 


1 

Z 


(F-4) 


Expression  F-4  can  then  be  used  to  approximate  the  far  field  potential 
flow  over  the  airfoil  of  chord  c  with  circulation  r  in  the  same  physical 
plane. 


Let  the  far  field  boundary  be  a  circle  of  radius  where 
Z  =  rf  e10,  0  <0  <  2-it .  Substitute  for  Z  in  Equation  F-4  and  obtain 


W,  (e)  =  e 


-la 


ela 
1 6rr~ 


-2i0 


CL  I  -i© 

4ir  r^ 


(F-5) 


u-iv , 

(F-6a) 
(F-6b) 

In  both  Equations  F-6,  the  first  term  is  the  freestream  velocity  contribution, 
the  second  term  is  the  doublet  contribution,  and  the  third  term  is  the 
vortex  contribution.  As  r^  approaches  infinity,  the  infinity  freestream 


Combine  the  exponential  terms  in  F-5  and  use  the  fact  that  U  = 
where  u  and  v  are  the  x  and  y  velocity  components,  to  obtain 

1  cl 

uf(o)  =  cos  a  -  2  cos  {a  -  2© )  +  sin  0 

1  Cl 

vf(0)  =  sin  a  +  fg^T2  sin  (a  -  20)  -  ^  cos  0 
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boundary  conditions  for  u  and  v  are  recovered.  Also,  the  doublet  term 
is  of  higher  order  in  1/r^  than  the  vortex  term  for  the  case  of  a  lifting 
airfoil  and  may  be  neglected. 


The  far  field  pressure  is  obtained  by  applying  Bernoulli's 
equation  for  irrotational  incompressible  flow  and  using  the  nondimensional 
form  for  pressure.  The  result  for  the  far  field  boundary  is 

Pf  (0)  =  1/2  {l-(uf2  +  vf2)}  (F-7) 

The  far  field  boundary  conditions  of  Equations  F-6  and  F-7  modify  the 
infinity  freestream  conditions  by  incorporating  effects  on  the  flow  of 
body  thickness  and  circulation  at  a  large  but  finite  distance  from 
the  airfoil . 

Transonic  small  disturbance  theory  for  slender  bodies  and  airfoils 
has  similar  expressions  for  the  far  field  conditions.  Klunker  (Reference 
101)  has  obtained  an  asymptotic  solution  applicable  at  large  distances 
from  thin  airfoils.  This  form  has  been  used  as  a  numerical  outer  boundary 
condition.  For  the  limiting  case  of  incompressible  two-dimensional  flow, 
the  doublet  and  vortex  strengths  are  compared  below  with  the  general  forms 
found  in  Equation  F-l . 

The  nondimensional  velocity  potential  doublet  in  the  far  field 
from  Klunker  becomes 

♦d  ’  s  lX F(s)ds  |F'8) 

where  F(s)  is  the  airfoil  thickness  function.  The  doublet  strength  is  thus 

2  yF(s)ds.  The  NACA  four  digit  airfoil  thickness  function  (Reference 

67)  is 

F(s)  =  (.2969  (£)*  -  .126  f  -  .3516(f)2  +  .2843  (f)3  -  .1015  (fM 

( F-9) 

where  c  Is  the  airfoil  chord,  f  is  the  maximum  thickness  as  a  fraction 
of  chord  and  s  is  the  chord  distance,  0  <_  s  <  c. 
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Substitute  the  thickness  function  Equation  F-9  into  the  doublet  strength 
expression  and  obtain  the  following  doublet  strength  from  transonic  small 
disturbance  theory: 

Uq  =  0.685  fc2  (F-10) 

For  the  case  of  a  NACA  0012  airfoil,  =  0.0822c2.  The  corresponding 

doublet  strength  in  the  simple  potential  flow  model  from  Equation  F-l  is 

2  2 

2 it  (^-)  or  .3972c  .  Klunker  points  out  that  the  doublet  may  be  neglected 

in  the  far  field  for  finite  circulation  flows  because  the  vortex  contri- 
1  1  2 

bution  is  of  order  (-  )  compared  with  an  order  (-  )  for  the  doublet. 
rf  rf 

The  vortex  potential  term  of  Klunker  for  the  case  of  incompressible 
two-dimensional  flow  is 

*V  =  i  {i  s9"  M  +  tan  '1^»  (F-m 

where  r  is  the  clockwise  circulation,  the  inverse  tangent  function  is 

defined  on  the  interval  (-  ^  j),  and  the  angle  within  braces  is  defined 

clockwise  from  the  negative  x  axis.  This  potential  function  is  identical 
with  the  vortex  function  in  Equation  F-l  where  the  clockwise  strength  is 
again  r  and  the  angular  orientation  is  the  right  hand  polar  coordinate  o 
defined  counterclockwise  from  the  positive  x  axis. 

Thus  the  far  field  potential  for  each  method  differs  only  in  the 
thickness  or  doublet  term.  Furthermore,  for  flows  with  finite  circulation 
this  difference  is  negligible  because  the  vortex  term  dominates. 
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APPENDIX  G 

NUMERICAL  INVISCID  FLOW  RESULTS 

A  numerical  solution  for  the  inviscid  flow  over  the  NACA  0012 
airfoil  is  obtained  for  comparison  with  both  the  Navier-Stokes  computation 
and  the  experimental  data.  The  numerical  algorithm,  developed  by  Thames 
(Reference  44),  which  solves  Laplace's  equation  for  the  stream  function 
in  two-dimensional  flow  is  used.  The  method  uses  SOR  iteration  to  solve 
the  system  of  equations  written  in  terms  of  body-fitted  coordinates. 

The  stream  function  value  on  the  body  can  either  be  specified  or  computed 
by  imposing  a  given  circulation  or  a  Kutta  condition. 

The  stream  function  equation  for  two-dimensional  flow  is 

i|i  +  ij;  =0  ( G-l  ) 

yxx  Tyy 

Rewrite  Equation  G-l  in  terms  of  the  transformed  variables  in  the 
computational  plane  using  expressions  in  Appendix  A  to  obtain 

a  +  ^nn  +  <*n  +  '  \  =  0  (G-2) 

The  boundary  condition  for  the  stream  function  on  the  outer  boundary,  which 
is  a  radius  of  ten  chords,  is  the  freestream  value 

<1^  =  yf  cos  <*L  -  xf  sin  aL  (G-3a) 

where  (x^,  yf)  are  coordinates  on  the  outer  computational  boundary  and 
is  the  angle  of  attack.  Hodge  and  Ghia  (Reference  82)  have  shown  by  an 

inviscid  analysis  that  this  boundary  condition  induces  an  error  in  lift 
of  less  than  1%  for  angles  of  attack  less  than  10°.  The  boundary 
condition  on  the  body  is 

=  constant  (G-3b) 

The  constant  can  be  specified  or  determined  by  imposing  a  given  circula- 
i  tion  or  a  Kutta  condition.  The  Kutta  condition  which  matches  the  upper 
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and  lower  trailing  edge  surface  velocities  by  extrapolation  is  selected. 

The  details  for  each  option  are  given  by  Thames  (Reference  44). 

The  derivatives  in  the  transformed  Equation  6-2  are  approximated 
by  central  difference  formulas  found  in  Appendix  B.  The  resulting 
system  of  equations  is  solved  by  SOR  as  described  by  Thames  (Reference 
44).  Solutions  are  obtained  for  0°,  5°,  7.5°,  9.5°,  and  11.5°  for  the 
NACA  0012  airfoil  section.  The  body  surface  pressure  coefficients  can  be 
calculated  from  the  stream  function  solution  as  follows.  The  body  pressure 
coefficient  is  defined  as 

Cp  *  1  -  (u2  +  v2)  (G-4) 

where  at  a  location  on  the  body  u  and  v  are  the  x  and  y  components  of  the 
flow  velocity  nondimensional ized  by  the  freestream  velocity.  Use  the 
definition  of  the  stream  function  to  obtain 

Cp  =  1  -  ((*y)2  +  (*/>  (G-5) 

Then  rewrite  the  Equation  G-5  in  terms  of  the  transformed  variables  and 
note  that  $  -  0  on  the  body  to  give 

Cp  =  1 '  ^  (V  (G’6) 

where  all  quantities  are  evaluated  on  the  body  surface.  The  body  pressure 
coefficients,  calculated  using  Equation  G-6  with  second  order  one-sided 
differences,  are  plotted  in  Figures  16,  17,  18,  19,  and  33.  A  typical 
streamline  plot  for  a  =  11.5°  is  shown  in  Figure  40.  The  numerical  in- 
viscid  solution  is  also  compared  with  an  analytical  solution  (Reference 
67)  which  uses  Theod orsen's  (Reference  12)  method.  The  excellent  agree¬ 
ment  is  seen  in  Figure  33  where  both  solutions  are  plotted. 
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